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SUMMARY 

A numerical algorithm and computer program are presented for solving the laminar, 
transitional, or turbulent two-dimensional or axisymmetric compressible boundary-layer 
equations for perfect-gas flows. The governing equations are solved by an iterative 
three-point implicit finite-difference procedure. The software, program VGBLP, is a 
modification of the approach presented in NASA TR R-368 and NASA TM X-2458, 
respectively. The major modifications are: (1) replacement of the fourth-order 

Runge-Kutta integration technique with a finite— difference procedure for numerically 
solving the equations required to initiate the parabolic marching procedure; 

(2) introduction of the Blottner variable-grid scheme; (3) implementation of an 
iteration scheme allowing the coupled system of equations to be converged to a speci- 
fied accuracy level; and (4) inclusion of an iteration scheme for variable-entropy 
calculations. These modifications to the approach presented in NASA TR R-368 and 
NASA TM X-2458 yield a software package with high computational efficiency and flexi- 
bility. Turbulence-closure options include either two-layer eddy-viscosity or 
mixing-length models. Eddy conductivity is modeled as a function of eddy viscosity 
through a static turbulent Prandtl number formulation. Several options are provided 
for specifying the static turbulent Prandtl number. The transitional boundary layer 
is treated through a streamwise intermittency function which modifies the turbulence- 
closure model. This model is based on the probability distribution of turbulent 
spots and ranges from zero to unity for laminar and turbulent flow, respectively. 
Several test cases are presented as guides for potential users of the software. 


INTRODUCTION 

A number of finite-difference and integral methods are currently available for 
numerically solving the two-dimensional, or axisymmetric, compressible boundary- layer 
equations. No attempt is made in the present paper to present a literature review of 
either solution techniques (ref. 1) or turbulence closure (ref. 2). Reference 2 
includes a tabular summary of 34 additional two-dimensional programs available as of 
1975. The purpose of the present paper is to present modifications of the algorithm 
and software presented in references 3 and 4 that render the approach more accurate, 
more efficient, and easier to implement. 

In the present approach, a coupled, iterative implicit finite— difference proce- 
dure, similar in many respects to that presented in references 5 and 6 for laminar 
flows, is used to solve the system of equations for laminar, transitional, or turbu- 
lent boundary-layer flows. The major modifications presented in the present approach 
as compared with references 3 and 4 are as follows: (1) replacement of the fourth- 

order Runge-Kutta integration technique used in reference 4 with a finite-difference 
procedure for numerically solving the equations required to initiate the parabolic 
marching procedure; (2) introduction of the variable-grid scheme proposed by Blottner 
in reference 7; (3) implementation of an iteration scheme allowing the coupled system 

of equations to be converged to a specified accuracy level; and (4) implementation of 
an iteration scheme for variable-entropy calculations. For most applications, because 
of the quasilinearization technique, the iteration cycle for constant-entropy calcula- 
tions can be omitted if a sufficiently fine grid distribution is chosen for the coor- 
dinate normal to the wall boundary. (See ref. 8.) The present program can be easily 
applied to any attached, compressible, perfect-gas (two-dimensional or axisymmetric) , 
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boundary- layer flow. Transverse- curvature terms are retained in the system of equa- 
tions with the option of being neglected if the user so desires. The program is also 
structured to allow the user the option of obtaining locally similar solutions. 

Options are provided for either two-layer eddy-viscosity or mixing-length 
turbulence-closure models. No attempt has been made to generalize the closure models 
to empirically include the effects of streamwise pressure gradient, wall curvature, 
wall roughness, wall mass transfer, and low Reynolds number. (See ref. 2.) The 
models are structured in subroutine TURBLNT such that the user can easily modify the 
formulation to best represent the specific type of flow under investigation. The 
static turbulent Prandtl number can be specified in one of three ways: (1) as a 

constant; (2) as an analytic function of the coordinate normal to the wall boundary; 
or (3) as tabular input from experimental data. 

The transitional region of the boundary layer is modeled through the streamwise 
intermittency function (ref. 9), which modifies the turbulence-closure model. 

Boundary- layer transition location and the extent of the transition (length of transi- 
tion region) can either be specified from experimental data or computed from empirical 
correlation equations. The laminar boundary- layer equations are recovered by equating 
the streamwise intermittency function to zero. 

Five test cases are presented in the present paper. These cases cover external 
and internal flows, including flows with wall mass transfer, transverse-curvature, 
and variable-entropy effects. The cases are designed to serve as guides for assisting 
potential users as they become familiar with program VGBLP prior to applying the pro- 
gram to their specific problems. 


SYMBOLS 

function , 26V/u^ 

constant, Au^/V 

coefficients in difference equation (45a) and defined by 
equations (B3) 

coefficients in difference equation (45b) and defined by 
equations (B4) 

a speed of sound 

a^,a 2 /a^,a^ coefficients in molecular viscosity relations (see eqs . (7)) 

Cf skin-friction coefficient 

Cn specific heat at constant presg^e. 

D damping term (eq. (17a)) 

F velocity ratio, u/u^ 


damping 


damping 


Al ,B1 ,C1 ,D1 
n' n' n n 




A2^,B2„,C2^,D2^ 




1 ir 


2 



heat-transfer coefficient 


h 
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k 

kn 


I 

I 

M 


N 




"Pr 


N, 


Pr , t 


^Re 


N. 


Re , r 


N. 


Re , s 
No 

Re,St,i 

%e,6* 

’^Re,0 


N, 


Re ,00 


N, 


St 


index used in grid-point notation (see eq. (41) ) 
geometric progression constant, Ar|j^_,_j^/ATi^ 
thermal conductivity 
eddy conductivity (see eq. (2c)) 

constant in eddy-viscosity model (see eq. (16a) ) 
constant in eddy-viscosity model (see eq. (16b)) 
constant in intermi ttency function (see eq. (17c)) 
constant in intermi ttency function (see eq. (17c)) 
constant in mixing-length model (see eq. (20b) ) 
reference length 
defined in equation (30) 
mixing length (see eq. (20a) ) 

Mach number 

grid-point index in S-direction 

number of grid points normal to wall boundary 

reference number of grid points normal to wall boundary (see eq. (42)) 
Prandtl number, Cpjj/k^ 

static turbulent Prandtl number (see eqs . (3)) 

unit Reynolds number, u^/v^ 

reference Reynolds number, 

Reynolds number based on s, ^e^/^e 

Reynolds number at transition, us. • /V 
^ e t,i' e 

Reynolds number based on displacement thickness, u^6*/v^ 

Reynolds number based on momentum thickness, u^0/v^ 
free-stream unit Reynolds number, 

Stanton number, h/ (c^pu) 
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n grid-point index in Y-direction (see fig. 1) 

P (1 ) f P (2 ^ { 3) defined in equations (50) 
p pressure 


Q 


( 1 ) 


.Q 


q 


R,Z 




r 


r 


o 


r 


s 


S,Y 


s 


^t,f 


s. . 
t, 1 

T 

t 


u 


V 


V 


V 


V 




Xi,X2,. 


Yi,Y2,. 


y 


(2)^g(3) defined in equations (50) 
heat- transfer rate 

body axis system with origin at stagnation point, where Z is positive 
downstream and R is positive radially outward (see fig. 2) 

gas constant (see eq. (6) ) 

radial body coordinate measured normal to Z-axis (see fig. 2) 
body radius (see fig. 2) 

radial coordinate of shock wave (see fig. 2) 

orthogonal boundary- layer coordinate system with origin at stagnation 
point, where S lies along the body surface and is positive downstream 
and Y is normal to the body surface and positive outward (see fig. 2) 

boundary- layer coordinate along S-axis (see fig. 2) 

end of transition (see fig. 2) 

beginning of transition (see fig. 2) 

static temperature 

transverse-curvature term (see eqs . (23)) 

velocity component in S-direction (fig. 2) 
friction velocity, ^T^/p 

transformed normal-velocity component (see eq. (26)) 
velocity component in Y-direction . 

_ . p * V ’ 

velocity component, v + 

P 

velocity component, ^ r 

..,X^ functions of grid-point distribution (see eqs. (A4) to (AS)) 

..,Y^ functions of grid-point distributions (see eqs. (A12) to (A17) ) 

boundary-layer coordinate along Y-axis (see fig. 2) 


y 


stretched y-coordinate (see eq. (15)) 


z 

a 


3 


r 

Y 

Y 


As, Ay 
As^ 
A?, An 

A* 

6 

6 * 

5* 

inc 


e 

e 

£ 


0 

0 

0, 

A 


U 

V 




P 

T 


match point for two-layer eddy-viscosity model 
axial body coordinate (see fig. 2) 
defined in equation (30) 
defined in equation (30) 

streamwise intermittency distribution (see eq. (38)) 
ratio of specific heats 

transverse intermittency distribution (see eq. (17c)) 
grid-point spacing, physical plane 
transition extent, ^t f ~ ^t i 

grid-point spacing, transformed plane (see fig. 1) 
defined in equation (50g) 
boundary- layer thickness 
displacement thickness 
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incompressible displacement thickness. 


eddy viscosity. 


^ 9u/9y 



(1 - F) dy 


defined in equation (5a) 
defined in equation (5b) 
static- temperature ratio, 
momentum thickness 
shock-wave angle (see fig. 2) 
defined in equation (40) 
molecular viscosity 
kinematic viscosity, y/p 

transformed boundary- layer coordinates (see fig. 1 and eqs . (22)) 

defined in equation (39) 
density 
shear stress 
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local surface angle (see fig. 2) 

X vorticity Reynolds number, ^ 

X maximum value of 

max '"^'m+1 

ip stream function (see fig. 2) 

Subscripts : 

aw adiabatic wall 

e based on boundary-layer edge conditions 

i inner region of turbulent layer 

m mesh point in ^-direction (see fig. 1) 

max maximum value 

n mesh point in g-direction (see fig. 1) 

o outer region of turbulent layer 

r reference quantity 

s shock 

t total condition 

w wall value 

free stream 
Superscript: 

j flow index; j = 0 for planar flow, j = 1 for axisymmetric flow 

An asterisk ( )* on a symbol denotes a dimensional quantity. 

A prime on a symbol denotes a fluctuating component. 

A bar over a symbol denotes the time average value. 

A coordinate used as a subscript denotes the partial differential with respect 
to the coordinate. (See eqs . (Al).) 
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PROBLEM DESCRIPTION 

This section presents the governing equations for compressible laminar, transi- 
tional, or turbulent boundary-layer flows together with the required boundary condi- 
tions, It should be noted that the system of equations can be found in numerous 
references (e.g., see refs. 2 and 3); however, for completeness, the equation set is 
presented in order to allow the user to modify the software if required. The alge- 
braic turbulence closure, transition location and extent, and transitional-f low- 
structure models are presented and briefly discussed; however, the reader interested 
in a detailed discussion of these models is referred to references 2, 3, and 8. 


Basic Partial Differential Equations 

Dimensional governing equations .- The mean turbulent boundary-layer equations can 
be written as follows: 

Continuity 


3s* 


(r*^p*u*) 


3y^ 


r*^p*\v* + 


p*v^ 


- 0 


(la) 


Momentum 


3u* 

3s* 


+ 1 V* + 


p*v* 3u* 
p* /3y* 


dp* 1 3 

ds* 3y* 


/ . 3u* ’ ’ 

V* ^ P*u*v* 


(lb) 


Energy 


+ 1-* + 


p*v* \ 3 


. ^(c*T*) 

dy* p 


. dp* 1 3 

“ ds* *j 3y< 


r*^ — ^A^(c*T*) 

c* dy* p 


3u* 1 


*j 9y* 


I I 

^-c*p*v*T* 


- p*u*v* 


3u* 

3y* 


(Ic) 


The mean turbulent equations are identical to those for laminar flow with the excep- 
tion of the correlations of turbulent fluctuating quantities. These correlations must 
be related to the mean flow in order to obtain a closed system of equations. In the 

I r i I 

present analysis, the apparent mass-flux term p*v* , the apparent stress term p*u*v* 

f I 

(Reynolds stress), and the apparent heat-flux term c*p*v*T* are modeled or repre- 
sented by a new velocity component v*, an eddy viscosity E* , and an eddy conduc- 
tivity k*, respectively, as follows: 
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V* 


e* 


k 


* 

T 


I I 


V* + 


-p* 



p*v* 

p* 


u*v* 

9u*/8y* 


* v*T* 

9T*/3y* 


(2a) 


(2b) 


(2c) 


The static turbulent Prandtl number is defined as follows: 


N 


Pr,t 


I I 

U* V* 


I f 

v*T* 


/8T*/3y* 

\8u*/3y* 


(3a) 


Equation (3a) can then be rewritten in terms of equations (2b) and (2c) as 


N 


Pr, t 


c*e* 

P 


(3b) 


In terms of equations (2) and (3), the governing differential equations can be 
written as follows: 

Continuity 


3 

3s* 


(r* 


j 


p*u*) 


+ 



p*v*) 


0 


(4a) 


Momentum 




. , . 3u* - 3u* 

P* lu* r + V* 


3s* 


3y^ 


ds* ^ j 3y* \ 3y* / 


(4b) 


Energy 




p*|u* (c*T*) + V* . (c*T* 


3s* p 


3y*"-p‘ >] - “* as* * ^‘(sy*) 


(4c) 
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e* appearing in equations (4) are defined as follows: 


G* 




(5a) 


e* 



N. 


Pr 


y* N 


Pr, t 


(5b) 


The function F appearing in equations (5) represents the streamwise intermittency 
distribution for the transitional-flow region. The F distribution assumes a value 
of zero for laminar flows, a value of unity for fully turbulent flows, and a range 
of 0 to 1 for the transitional region of flow. The variation of F in the transi= 
tional region depends upon the statistical growth and distribution of turbulent 
spots. The model used to represent F is discussed in a subsequent section of the 
present paper. 

The system of equations is closed by the addition of the perfect-gas law and a 
viscosi ty-temperature relation. The perfect-gas law is written as 


p* = p*R*T* 

g 


( 6 ) 


Two viscosity- temperature relations are provided: (1) the Sutherland law 




a* (T*) 


T* + a 


* 

2 


(7a) 


and (2) the power law 


y* = a* (T*)^"^ 


(7b) 


The pressure-gradient term in equations (4) is replaced by the Bernoulli 
relation 


dp* 

ds* 


~p*u* 

e e 


du" 


ds* 


( 8 ) 


for constant entropy flows; however, for variable entropy flows the value of dp*/ds* 
is explicitly retained in the equation system. 

The equations are rewritten in nondimensional form where the nondimensional 
variables are as follows: 
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u = u*/u* 
r 


V = v*/u* 
r 


P = P*/(p*u^*) 


p = p*/p* 


»p = fp* /T* 

r 


y = y*/y* 






( 9 ) 


where all dimensional lengths are nondimensionalized by a reference length L* . The 
reference values of density and velocity are taken to be those of the free stream, the 
reference temperature is taken to be reference viscosity is the 

viscosity obtained from either equation (7a) or (7b) evaluated at the reference 
temperature . 

Nondimensional governing equations .- The nondimens ional equations are as 
follows : 

Continuity 


-^(r^pu) + — (r^pv'*') = 0 


9s 


3y 


( 10 ) 


Momentum 




9u + 9u 
p u ^ + v+ — 

9s 9y 


ds 95 V By/ 


( 11 ) 


Energy 




9t 

9y 


ds 


1 9 / 9t\ ~/9u\^ 

— — e — 

9y \ 9y/ \9y/ 


( 12 ) 


Equation of State 


P = 



(13) 
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QUALITY 


or 


= a.3/2/ 


1 + a. 




y = T 




+ a. 


y = T 


( 14 ) 


where 

y = ^)I^Re,r 

= ^{^Re,r > 


Algebraic models are used to close the system of equations, 
(1) a two- layer eddy-viscosity model (KODVIS = 2 ) , and 
(KODVIS - 1). 

Two- layer model 




Turbulence closure . - 
Two options are provided: 
(2) a mixing-length model 


The equations describing the two-layer model are as follows (see ref. 8) : 



D* 9 

^(kiy*D)2 


9y^ 


)*Y 

y* 2 e inc 


where 


D = 1 - exp(-y*/A*) 



(0 < y* < y*) (16a) 


(y* < y*) (16b) 


(17a) 


(17b) 


and 
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Y = 



(17c) 


The boundary layer thickness 6 in equation (17c) is defined as the distance normal 
to the wall boundary where = 0.995. The empirical constants to k^ are 

assigned values of 0.4, 0.0168, 5.0, and 0,78, respectively. 

The location of the boundary separating the two layers y* is defined from the 
continuity of eddy viscosity; that is, where ^ 



Mixing-length model 

A mixing— length option (KODVIS = 1) is provided for those interested in utilizing 
experimental mixing-length data. The eddy-viscosity distribution across the boundary 
layer can be written as follows: 


_ P* ^*2 3u* 
y M* 9y* 


(19) 


where the mixing length £* can be expressed as 



(20a) 


An analytic formulation is provided in subroutine TURBLNT for f [ x J as follows (see 
ref. 10) : / 


f 



= kc 


tan h 



(20b) 


where k^ is assigned 
values of k^, k^, ... 
to program VGBLP. 


the value of 0.108. It should be noted that the assigned 
^5 definition of 5 can be modified through input 
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The eddy conductivity (eq. (2c)) is modeled as a static turbulent Prandtl number 
(eq. (3a) ), Three options are provided in subroutines TURBLNT for ^pr,t' ^ 

constant, for example Np^ ^ = 0.95 (KODPRT = 1); (2) the Rotta (ref. 11) 
distribution (KODPRT = 2) 


N 


»Pr,t>w 


Pr, t 



(21) 


and (3) a distribution Np^ ^ specified in tabular form from experimental data 

(KODPRT = 3) - 


Transformed plane .- The system of governing equations is singular at s = 0. 

The Probstein-Elliott (ref. 12) and Levy-Lees (ref, 13) transformation is used to 
remove this singularity as well as to reduce the growth of the boundary layer as the 
solution proceeds downstream. This transformation can be written as follows: 


K (s) 



2 J 

p u y r 
e e e o 


ds 


(22a) 


n (s,y) 


pur 
e e o 




(22b) 


where the parameter t appearing in equation (22b) is the transverse-curvature term, 
defined as 


t 



(23a) 


or, in terms of the y-coordinate , as 


t = 1 + cos cb (23b) 

r 

o 

The relation between derivatives in the stretched physical (s,y) and transformed (C^H) 
coordinate system is as follows: 


fd \ _ 2j/3 \ ^ 

\8sj~ ^e^e^e^o \9sj\3ri/^ 

3 \ _ ^e%^o^ / p V 3 \ 

dyJs ~ M 


(24a) 


(24b) 


Two parameters F 
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and 0 are introduced and defined as 


F 


0 



T 

T 

e J 


( 25 ) 


together with a transformed normal velocity 


V 


2? 


2 j 

p u y r ^ 
e e e o 



pvr^ t^ 
o 

+ 

\[^ 


(26) 


The governing equations in the transformed plane can then be written as follows: 
Continui ty 


3v 9 f 

Tp- + 2^ ^ + F = 

9n ^ 


(27) 


Momentum 


^ „ II _ ^L2jy- ^'] ^ 


2^F 9^ + V g^lt 7e ) + B(F 0) 0 


Energy 


prp M + V 3 


f) - = 0 


where 


(28) 


(29) 



(30) 
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By using the viscosity relations (eqs. (14)) and the equation of state (eq. (13)), the 
parameter t can be written as follows: 

Sutherland law 


I - •So 


* V^e \ 
0 * a^/Tj 


Power law 


I = 


(0) 


^4 


-1 


Ola) 


(31b) 


The transverse-curvature term can be written in terms of the transformed vari- 
ables as 


t 



2\[2^ cos (|) 

p u r^^ /nT 
e e o V Re,r 



0 



(32) 


The physical coordinate normal to the wall is obtained from the inverse transforma- 
tion; namely. 


y 



cos (j) 


-1 



-f* 


2^|^^ cos 4. 

I 0 dn 

p u r Jo / 

"^e e o V Re,r / 


(33) 


The boundary conditions in the transformed plane are as follows: 
Wall boundary 


F(C,0) =0 1 


V(^,0) = V„(C) 


or 



(34a) 
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Edge boundary 


F(^.ne) = 1 I 

0(C,ne) = 1 


(34b) 


The boundary condition at the wall for the transformed V component can be related to 
the physical plane as (see eq. (26)) 


Vw = 




y r 
e o 


3 \P^^, 


e e> 


(35) 


Transition 

Transition location Because of the large parameter space influencing transition 
to turbulence (refs. 14 to 17)^ it is not possible to predict with assurance the loca- 
tion of transition for general flows. However, for certain classes of geometry 
(e.g., flat plate, cone, etc.), empirical correlations are available. These empirical 
correlations can be used with confidence provided one realizes that a probable range 
of transition locations is being predicted. In program VGBLP either the transition 
location (SST) or the stability index at transition (SMXTR; see ref. 8) must be 
specified; however, any correlation can be directly incorporated into the program. 

Transition extent .- The assumption of a universal intermi ttency distribution 
implies that the transition- zone length (transition extent) can be expressed as a 
function of the Reynolds number at transition reference 9 it is shown, 

for the transition data considered, that the data are represented on the average by 
the equation 


%e,As^ 


5/n 


Re,s^^i 



(36) 


0 

where ^Re,As ~ ~^^t f ” ^t i^ ’ location of the end of transition £ can 

t 0 ^ ^ ^ 

then be obtained directly from equation (36) as follows: 


"t, f 


^t,i 


+ 5N 


_-l 

Re 




(37) 


where is the local unit Reynolds number, 

In program VGBLP the extent of the transition region (s^. ^ ■) can be 

specified in one of two ways: (1) from equation (37) (KTCOD =1); or (2) from the 

specification of ^t,f^^t,i ot>tained from experimental data (KTCOD =2). It should 
be noted that the program can be easily modified to include any desired correlation 
or equation in place of equation (37) . 
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Transition- flow structure . “ The parameter T (eqs. (5)) is the streamwise inter- 
mittency function which models the turbulent spot distribution in the transitional 
region. The parameter is a function of the s-coordinate only and is defined (see 
ref. 9) as follows: 



r(^) = 1 - exp (-0. 412^^) 


(38) 


where 





and 


(39) 


A = 




(40) 


It should be noted that F - 0 for laminar flow, F = 1 for fully turbulent 
flow, and F ranges from 0 to 1 for the transitional- flow region. Equations (27) 
to (29) reduce to the classical laminar boundary-layer equations when F is set 
to zero . 


SOLUTION TECHNIQUE 

The system of governing equations (eqs. (27) to (29)) is parabolic and, therefore, 
can be numerically integrated by a marching procedure in the streamwise direction. In 
order to cast the equations into a form in which th^' marching procedure can be imple- 
mented, the derivatives with respect to ^ and r\ are replaced by finite-difference 
quotients. The method of linearization and solution used in the analysis closely 
parallels that of references 5 and 6. 


Finite-Difference Mesh 

Because of the magnitude and variation of the gradients of the dependent variables 
(3F/3y; 30/3y) near the wall boundary for turbulent flow, it is computationally 

inefficient to use equally spaced mesh points in the y-coordinate . This problem can 
be alleviated by selecting a variable mesh-point distribution such that An > 1 
where the distribution in the wall region is chosen sufficiently small to resolve all 
gradients. One approach to grid specification is to use a geometric progression 


An^ = (k)^ ^ Ai1]_ (i = 2, 3, 4, N) (41) 

where k is defined as the geometric progression constant An . A schematic of 
such a grid is presented (not to scale) in figure 1. 

Blottner (ref- 7) introduced a variable-grid scheme that is more computationally 
efficient than the differencing scheme and mesh distribution used in references 3 
and 4. The Blottner variable-grid scheme (ref. 7) and the Cebeci-Keller box scheme 
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(ref. 18) appear to be two of the most promising schemes currently available in the 
literature for solving the boundary- layer equations. The Cebeci-Keller box scheme, 
although efficient for two-dimensional flow, yields block- tridiagonal matrices and 
requires greater computational effort than the simple-tridiagonal matrices of the 
Blottner variable-grid scheme. Blottner has shown (ref. 7) that the variable-grid 
scheme is as accurate as the box scheme for the two-dimensional boundary- layer equa- 
tions and, furthermore, that large values of the geometric progression constant can be 
used for the normal mesh- point distribution, provided the variable grid is interpreted 
as a coordinate transformation. In reference 8 it was shown that k values on the 
order of 1.04 could be used for turbulent flows. In reference 19 it was shown that an 
accuracy requirement on i^ of 1 percent required approximately 220 mesh points 
normal to the wall with k - 1 . 04 . In order to increase the computational efficiency 
of such schemes one can reduce the number of mesh points while simultaneously increas- 
ing the value of k; however, this approach generally results in unacceptable levels 
of accuracy. Blottner (ref. 7) demonstrated that with the variable-grid scheme satis- 
factory results could be obtained with approximately 20 mesh points for k = 1.82. 
Vatsa and Goglia (ref. 20), using the method of reference 4, showed that the variable- 
grid scheme proposed by Blottner (ref. 7) could reduce the number of grid points from 
approximately 201 to 61 for a specified 1-percent accuracy in wall shear stress. They 
also showed that for most applications one could obtain reasonably accurate solutions 
for turbulent boundary layers with as few as 25 to 30 mesh points, as compared to 
approximately 200 mesh points for the method of reference 3. 

Blottner (ref. 7) introduced a modified definition for the geometric progression 
constant such that 


k 


N-1 



(42) 


where k is now defined as the conventional geometric progression constant for 
mesh points normal to the wall. Using equation (41), one obtains 





(43) 


which when combined with equation (42) yields the following for the p-mesh 
distribution : 


= "lel 


_ (i-1) (N -1)/(N-1) 

ik) ^ ^ 

N-1 

(k) ^ - 1 


(44) 


Two options are provided in program VGBLP: (1) specify Uitiax' ^ 

(IGEOM = 1); or (2) specify Nf and Ap. (IGEOM =2). Of these two options, 

max j. 

it is recommended that the user select the first option (IGEOM - 1) where the value 
of k is computed from equation (42) for specified values of k and N^. Typical 
values for k and are 1.5 and 25, respectively, for N ^ 41. (See ref. 20.) 

It is obvious that the larger the selected value of N and the smaller the value 
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of k - 1 (k ^ 1), the more accurate the solution. Since program VGBLP is very effi- 
cient in terms of computer resources, it is suggested that potential users of the 
program perform a series of numerical experiments over a range of k and N values 
in order to gain experience with the procedure. If the second option (IGEOM = 2) is 
selected (specify the user is cautioned to exercise care in 

selecting the number of grid points. If transitional or turbulent flow occurs in a 
given problem, the laminar region of the boundary layer is calculated with the value 
of k used for the turbulent region; that is, for a given solution, k is invariant. 


Difference Equations 

Three-point implicit difference relations (see appendices A and B) are used to 
reduce the transformed momentum and energy equations (eqs. (28) and (29), respectively) 
to finite-difference form. The difference quotients produce linear difference equa- 
tions when substituted into the momentum and energy equations provided truncation 
terms of the order A^^ and Ar)j^_jl^ Aq^^ are neglected. (It should be noted 

that the truncation term for d^F/dr]^ is of the order (Arij^«2 " An^^).) The resulting 
difference equations may be written as follows; 

Momentum equation 


®^n^m+l,n , n+1 ^ , n-1 




(45a) 


Energy equation 


^^n^m+l,n-l ®^n^m+l,n ''^n^m+l,n+l ^^n®m+l,n-l 




(45b) 


The coefficients Al^, Bl^, ..., Gl^ and A2^^, B2^, G2^ (see appendix B) are 

functions of known quantities at stations m and m-1. It is important to note that 
equations (45) are coupled through the dependent variables F and 0; however, the 
dependent variable V does not appear explicitly as an unknown at station m+1 . The 
variable V is uncoupled from the system because of the particular way that the non- 
, . 3 f ^ 30 

linear terms V ^ and V are linearized. (See eq. (A23).) 


Solution of Difference Equations 

The system of difference equations (eqs. (45)) represents a set of 2 (N - 1) 
linear algebraic equations for 2 (N - 1) unknowns. The boundary conditions to be 
used with the difference equations are specified in equations (34). The 2 (N - 1) 
linear algebraic equations may be written in tridiagonal matrix form; consequently, 
an efficient algorithm (Thomas algorithm) is available for simultaneous solution. 
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The simultaneous or coupled*-solution technique is presented in appendix B of 
reference 5; however, because of differences between the present work and that pre- 
sented in reference 5, the solution technique is discussed here in some detail. 

Because of the special form of equations (45) , the following relations exist 
(see ref. 21) : 




^ r.(2) 

^(3) 


^m+l,n-l 

- ^m+l,n-l 

+ P F 

m+1 , n-l^m+1 , n 

^m+l,n-l®m+l,n 

(46a) 

®iti+l,n-l 

(1) 

(2) 

(3) ^ 


ym+l,n-l 

2j^f^+l,n-l m+l,n 

®m+l ,n-l m+1 , n 

(46b) 


Next, equations (46) are substituted into equations (45) to obtain the following 
relations : 


* F + Fl *0 Cl * — Cl P 

m+l,n m+l,n m+l,n m+l,n m+l,n '^m+l,n m+l,n+l 


- FI 0 

m+1 , n m+1 , n+1 


T5T* p 4 , p9* 0 C9* — C9 F 

m+l,n m+l,n m+l,n m+l,n m+l,n "^m+l,n m+1, n+1 


— F9 0 

m+1 , n m+1 , n+1 


where 


and 


Bl 


m+1 , n 


Bl ,, + Al , + D1 

m+l,n m+l,n m+l,n-l m+1 , n-m+1 , n-1 


El 


m+1 , n 


El + Al n + D1 

m+l,n m+l,n m+l,n-l m+1 , n'^m+l , n-1 


G1 


m+1 , n 


- Al P<^^ - D1 

m+l,n m+l,n m+l,n-l m+1 , n^m+1 , n-1 


B2*_^- 
m+1 , n 

E2*.t ^ 

m+i , n 


( 2 ) ( 2 ) 

B2 + A2 , + D2 , 9 ^ 

m+l,n m+l,n m+l,n-l m+l,n m+l,n-l 


E2n,+l,n + ^2^+l,nPm+l,n-l + ^2^+1 , n^m+l , n-1 


'* = G2 - A2 _ Ts'p 

m+l,n m+l,n m+l,n m+l,n-l ^m+1 , n^m+1 , n-1 


(47a) 

(47b) 

(48a) 

(48b) 

(48c) 

(48d) 

(48e) 

(48f) 
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The unknown values of F and 0 at station m+l,n are obtained from equations 
as follows: 


F = F + P^^^ 0 

m+l,n m+l,n m+l,n m+l,n+l m+l,n m+l,n+l 


where 


0 = F + 0 

m-(-l,n ^m+l,n ^m+l,n m+l,n+l ^m+1 , n^m+1 , n+1 


^m+1 , n (^^m+1 ,n^^m+l , n ^^m+1 , n^^m+1 , n)^m+l , n 


= (e1*^, C2 - E2*^, Cl \A 

m+l,n \ m+l,n m+l,n m+l,n m+l,ny i 


^m+1 , n (^^m+1 , n^^m+1 , n ^^m+1 , n^^m+1 , n) 


m+1 , n 


= «+l,nGCl,n " , n^Cl , n)Cl , n 


o ~ ( on — c* y ^A^ 

^m+l,n ' m+l,n m+l,n m+l/n m+l,n/ m+l,n 


^m+1 , n (^^m+1 , n^^m+1 , n ^^m+1 , n^^m+1 , n) ^m+1 , n 


and 


m+l,n 


(®^m+l,n^^m+l,n ®^m+l , n^^m+1 , n) 


Next, equations (46) are rewritten as follows (where n = n + 1) 


F +P^2) p ^ p(3) Q 

m+l,n m+l,n m+l,n m+l,n+l m+l,n m+l,n+l 


0 = + n^2) p . „(3) „ 

m+l,n ^m+l,n ^m+1 , n^m+1 , n+1 ^m+1 , n m+1 , n+1 


The "no'-slip" boundary condition i ” applied at the wall boundary to 

obtain the values of i ^here i = 1, 2, 3; that is. 


p(l) =p(2) = p(3) ^0 

m+1,1 m+1,1 m+1,1 


( 47 ) 

(49a) 

(49b) 

(50a) 

(50b) 

(50c) 

(50d) 

(50e) 

(50f) 

(50g) 

(51a) 

(51b) 
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The thermal condition at the wall boundary can be specified in one of two ways: 

(1) specified wall-temperature distribution (KODWAL = 1) ; or (2) specified heat- 
transfer distribution (KODWAL = 2) . For a specified wall- temperature distribution it 
can be seen directly from equation (51b) that 


( 1 ) 

^m+1,1 


= 0 




m+1 , 1 


■^m+1 , 1 


(3) 

■^in+l , 1 


= 0 


(53) 


The case in which a heat-transfer distribution is specified presents a somewhat more 
difficult problem; however, this class of flows is often of interest; for example, 
adiabatic flows where q* = 0, 

The heat transfer at the wall boundary can be written in the transformed plane 
as follows (see ref. 3) : 


^m+1 , 1 


*^^*2 , 
^r r V Re,r 

'p u T ]A r^\ 
e e*^e o \ 

* \ 

0^ 


I 

™+l'l\3n/m+l,l 


(54) 


Then, for a specified value of th® gradient of 0 can be obtained directly 

as follows : * 


\9n/m+l,l ^m+ 1,1 * *2 


^r^r , r V’e’^e'^e^e^o/ ^ 


(i) 

\^/m+l,l 


(55) 


For the grid-point spacing used in program VGBLP, the gradient of 0 evaluated at 
the wall, by using a three-point relation, is as follows: 


^9n/m+l,l 


[l - (1 II - 

k(l + k)An, 


0 


m+1 , 3 


(56) 


Equations (55) and (56) then yield the following expression for 0jn+i i = 
k(l + DAn^/^v 


0 


^ 0 ., „ + — ^ 


0_ 


1 - (1 + k) ^ ''^'^^m+1,1 1 - (1 + k)^ m+1, 2 1 _ (1 + k)2 m+l,3 


(57) 


where / is evaluated from equation (55) . Equations (45) are next written 

at the m+1, 2 point to obtain two equations in terms of ^j^i+1 n ®m+l n 
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n = 1, 2, 3. (Note that ^ = 0.) The quantity ^ is then eliminated from 

these two equations to obtain one equation in terms of 2 ®m+l n 

n = 1, 2, 3. The quantity 3 next eliminated through use of equation (57) 

to obtain the relation 


0 = 5^^^ + 5^2) p ^ 5O) g 

'^m+1,1 ^m+l,! + ym+l,l^m+l, 2 ^m+l , l^m+1 , 2 


(58) 


where 


^m+1,1 


^(2) 

*^m+l , 1 




[(C2)1G1) - ^ [(C2KF1) - (c1)(F2)]^^j_2[’'“ * 


^m+l/2 


[(Cl) (B2) - (C2) (Bl)] . 


A 


m+1 , 2 


[cD (E2) - (C2)(E1)^^^^ 2 + [ci) (F2) - (C2) (Fl^ ^{1 + )c) 


m+ 1 , 1 


(59a) 


(59b) 


(59c) 


m+1 , 2 


and 


A 


m+1 , 2 


= |[C2)(D1) - (Cl) (D2)] + [C2) (FI) - (Cl) (F2)] [l - (1 + )<)^| 


(59d) 


m+1 , 2 


By comparing equations (51b) and (58) , it is observed that 


(i) 

"^m+l , 1 


= Q 


(i) 

m+1 , 1 


(i = 1, 2, 3) (60) 


which completes the desired boundary condition for the case of a specified heat- 
transfer distribution along the wall boundary. The temperature at the wall is 
obtained from equation (57) once ®rn+X 2 ^m+1 3 known. 

The quantities ®m+l n i = 1, 2, 3 (see eqs. (51)) must 

first be determined across the boundary layer at the m+1 station where 
n = 1, 2f . . , , N. These quantities are calculated by the following procedure: 

(1) Perform the following steps at the first grid point away from the wall 
(n = 2) : 


I II 
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(a) Calculate ^^2' equations (B3) 

(b) Calculate ^'^2* *•*' ^^2 ^'^om equations (B4) 

(c) Using the results from steps (a) and (b) and the boundary conditions 

(eqs. (52) and (53) or (59)), calculate Bl^, B22/ El*r E2*, Gl*, 
and G22 from equations (48) 

(d) Using the results from steps (a) to (c) , calculate and 

where i = 1, 2, 3 from equations (50) 


(2) The procedure outlined in step (1) is now repeated at the second grid point 
off the wall (n = 3) by using the results obtained at n = 2. This procedure is 
repeated until the entire boundary layer is traversed (n = N) and all values of 

^m+l,n ^m+l,n determined where i = 1, 2, 3 and n = 2, 3, 4, . . . , N. 


,(i) 


,(i) 


where i = 1, 2, 3 and 

where n = N-l, N-2, 


(3) Using the values of and O 

Tn+l,n ^m+l,n 

n=2, 3, 4, the values of F and 0 

m+l,n m+l,n 

are calculated from equations (49). It should be noted that F . , and 0 , 

m+l,N m+l,N 

specified edge boundary conditions (eqs. (34b)), The wall-boundary values of F 
and 0 are obtained from equations (34a), or from equation (57) for the case of , 
specified wall-boundary heat-transfer distribution. 


2 

are 


the values of F 


m+1 , n 


and 0, 


only to determine V 


m+1 , n 


m+1 ,n 
for all values of 


At this point in the procedure, 
are known for n = 1, 2, . . . , N, and it remains 

n to complete the first iteration. 


(4) The continuity equation (eq. (27)) is solved numerically for the N - 1 


unknown values of V 


m+1 ,n 


as follows: 


^m+l,n ^m+1,1 



dn 


(61) 


where ^ represents the wall-boundary condition V^. (See eq. (35).) The 

trapezoidal rule of integration is used to numerically solve equation (61). 

(5) The solution is now checked for convergence where the convergence criterion 
is as follows (q is iteration index) : 


DIF 


1 


OF/3n)^ 

OF/9n)^"^ 


m+1 , 1 


(62) 


If DIF ^ CONV, station m+1 is declared converged and m is incremented by 1. 
During the development of program VGBLP, a global convergence check was initially 
applied to each of the three dependent variables over all values of n. Global con- 
vergence was declared once all three variables were converged over all values of n; 
however, it was observed numerically that a local check on the gradient of F at the 
wall was a sufficiently accurate criterion. Consequently, the logic and storage 
requirements for the global check were removed from the software. 
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Initial Profiles 

A major change in the present approach compared with that of references 3 and 4 
is the technique for numerically generating the initial profiles at ^ = 0. These 
initial values are required to initiate the three-point implicit marching procedure. 

In references 3 and 4, the equations at the initial station = 0) were numerically 

solved using a fourth-order Runge-Kutta scheme for an equally spaced grid 

(An = n max^ - 1)) in the q-direction . The converged solution on the equally spaced 

grid was then interpolated onto the variably spaced grid = k . This pro- 

cedure introduced some interpolation error into the initial profiles which, although 
decaying in could cause oscillations in the neighborhood of the stagnation point 

for blunt-body flows. These oscillations, if they occurred, made it difficult to 
accurately determine parameters such as lim q^ and lim t^. Another problem some- 

s->0 s^O 

times encountered by users of the approach presented in reference 4 was the sensi- 
tivity of the convergence of the initial solution to the selected initial values 

required to initiate the Runge-Kutta integration. To ensure 

n=o 

convergence, the user would often have to make several trial runs before the initial 
guesses were sufficiently close to the converged values- 

In the present approach these two problems (interpolation and convergence) are 
completely eliminated. The locally similar form of the momentum and energy equations 
is numerically solved in finite-difference form for the variable grid used in the 
marching procedure. Initial value guesses and interpolation procedures are not 
required. The momentum and energy equations are of tridiagonal form and are easily 
solved. The continuity equation is uncoupled from the system and solved using the 
trapezoidal rule of integration. The initial profiles are generated in subroutine 
SIMILAR. 

The marching procedure requires evaluation of the ^-derivatives at two backward 
points; consequently, the first solution station SS(1) + SS{2) in the ^-coordinate 
requires special attention to assure local accuracy in the neighborhood of the stag- 
nation point. For flows with a stagnation point (IBODY = 1) , the profiles at 
C = are reflected about the stagnation point to impose symmetry. For flows 

without a stagnation point (IBODY =2), it is assumed that the solution at ^ = A^l 
is identical to that at ^ = 0 (s - 0) . As a result of this approach, the solutions 
are physically correct and merge smoothly with the downstream marching solution 
(^ > A^j^) . The quantity A^2. evaluated using Simpson’s rule for stagnation-point 

flows . 


3f 


and — 
3n 


Evaluation of Wall Derivatives 

The shear stress and heat transfer at the wall are directly proportional to the 
gradient of F and 0 evaluated at the wall, respectively. By using G to repre- 
sent a general quantity, where \ is the value of G at m+1 evaluated at the 

wall, the four-point difference scheme used to evaluate derivatives at the wall is 
given as , 


(i), 


m+1 , 1 




(63) 


:! II 


26 


OFPOOft^AuT? 

where the coefficients are defined by the following relations: 


and 


(1 + k + k^)^[k(l + k) - l3 + (1 + k) 
(1 + k) (1 + k + k^)k^ Arij^ 


(1 + k + k^) 

An^ 

(1 + k + k^) 
(1 + k)k^ Ari2_ 


Y 


10 


1 

(1 + k + k^)k^ Arij^ 


(64a) 


(64b) 


(64c) 


(64d) 


For the case of equally spaced grid points in the q-direction (k = 1) , equa- 
tion (63) reduces to the familiar four-point relation: 


/m+l , 1 


6 An^^^'^m+1,1 ^®^ra+l,2 ^^m+1,3 


^^m+1,4^ 


(65) 


PROGRAM DESCRIPTION 

Program VGBLP is run on the Control Data CYBER 170 series computers under the 
NOS 1.4 operating system at the Langley Research Center. 


Array Dimensions 

Program VGBLP and subroutines TABLE, VARENT, TURBLNT, MESH, SIMILAR, and SOLVE 
use the variable-dimension capability of the preprocessor at the Langley Research 
Center. This capability allows the user to designate the minimum storage require- 
ments for a given problem. If the preprocessor capability is not available at the 
user s installation, the dimension statements can be modified by inserting the 
required dimensions in place of their equivalent designations (UPDATE, MODIFY, etc.) 
in accordance with the following definitions for program VGBLP and its subroutines. 

Program VGBLP 

As^ sig ned va lue 

1 - (KODPRT ^ 3); NUMBl - (KODPRT = 3) 

IE + 2 


Variable dim e nsio n 
JI 
JK 
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1 - (constant entropy); lENDl - (variable entropy) 
INTEGER (s„^^/PRNTINC + IPRNT) 

ITlaX 

INTEGER (s^^ /PROINC + IPRO) 

IilclX 

TENDl + 1 


Siobroutine TABLE 


JH See definition in VGBLP 

JJ NUMBER 


Subroutine VARENT 


JL 


See definition in VGBLP 


Subroutine TURBLNT 


JI 

See 

definition 

in VGBLP 

JK 

See 

definition 

in VGBLP 


Subroutines MESH, SIMILAR, SOLVE 

JK See definition in VGBLP 


Input Description 

Standard CDC NAMELIST is used for all data input. Program VGBLP reads input 
under $NAM1 . Subroutine TABLE reads input under $NAM2. For cases where the variable 
entropy option is required (lENTRO = 2 ), subroutine VARENT reads input under $NAM3. 

Input/output flexibility is provided to the user wherein either the International 
System of Units (SI, KODUNIT = 1) or the U.S. Customary Units (U.S., KODUNIT = 0) can 
be used. The required input and resulting output data are listed in the following 
sections with appropriate dimensional units. The SI Units are listed first, followed 
by U.S. Units in parentheses. If no units are listed, the quantity is nondimens ional . 

Input data for $NAM1 . - 


Variable name 


CONV 

CONVE 

DETAl 


Variable description 

Convergence criterion for boundary- layer solution; 
DEFAULT = 1 X 10"^ 

Convergence criterion for variable entropy 
iteration; DEFAULT = 1 x 10"^^ 

(see fig. 1) 
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FT 

G 

GLAR 

IBODY 

IE 

lENDl 

lENTRO 

IGAS 

IGEOM 

IPRNT 

I PRO 

ITMAX 

lYINT 

J 

KODAMP 

KODE 

KODPRT 


1.0 - nonsimilar solution; 0.0 - locally similar 
solution; DEFAULT = 1.0 


Y; DEFAULT =1.4 


y/5 array corresponding to PRTAR array, used only 
if KODPRT = 3; NUMBl values 

1 - flows with stagnation point; 2 - flows without 
Stagnation point 

Number of mesh points in r}~coordinate 

Number of steps in S-direction 
lENDl 

^max " ^ SS(m) where SS(1) =0 

m=l 


1 constant entropy; 2 — variable entropy; 

DEFAULT = 1 

1 - Sutherland's viscosity (see eq. (7a)); 

2 - power-law viscosity (see eq. (7b)); 

DEFAULT = 1 

1 - specify XEND, IE, XK; 2 - specify XEND, IE, DETAl 

Number of specified wall-value printouts desired 
other than those determined by PRNTINC; 

DEFAULT = 0 

Number of specified profile printouts desired other 
than those determined by PROINC; DEFAULT = 0 

Maximum number of iteration cycles, 1 < m < lENDl , 
for variable-entropy calculations; DEFAULT = 3 

1 - normal intermittency function, y set to 1.0; 

2 - normal intermittency function, y calculated 
from equation (17c); DEFAULT = 1 

j; 0 - two-dimensional; 1 - axisymmetric 

1 - local values used in equation (5b) for damping; 

2 wall values used in equation (5b) for damping; 
DEFAULT =2 

1 - both laminar and turbulent profile prints are 
desired for diagnostic reasons once flow is 
turbulent; 0 - otherwise; DEFAULT = 0 

1 - constant 2 - Rotta distribution (see 

eq. (21)); 3 - tabular, Np^ . = f(y/6); 

DEFAULT = 1 / c 
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QUAUTY 

KODUNIT 


0 - all dimensional input and output in U.S. Units; 
1 - all dimensional input and output in SI Units; 
DEFAULT = 0 

KODVIS 


1 “ mixing-length model; 2 - two-layer eddy- 
viscosity model; DEFAULT = 2 

KODWAL 


1 - specified wall temperature distribution; 

2 - specified wall heat- transfer distribution; 
DEFAULT = 1 

KTCOD 


1 - transition extent calculated from equa- 
tion (37) ; 2 - transition extent specified as 
TLNGTH; DEFAULT = 2.0 

NAUXPRO 


1 - auxiliary profile prints are desired (see out- 
put description); 2 - otherwise; DEFAULT = 2.0 

NITMAX 


Maximum number of iterations allowed at any given 
station; DEFAULT = 1 

NUMBl 


Number of values read into PRTAR and GLAR arrays if 
KODPRT = 3 

PHII 


-l/ds\ . 

Opening angle of body at s - 0, tan u^g 

PR 


; DEFAULT = 0.72 
Pr 

PRNTINC 


Incremental s*-value for which wall-value printouts 
will be made, m (ft); DEFAULT = 0.1 

PRNTVAL 


Array of IPRNT specified s*-values for which wall- 
value printouts are desired, m (ft) 

PROINC 


Incremental s*-value for which profile printouts 
will be made, m (ft); DEFAULT = 1.0 

PROVAL 


Array of IPRO specified s*-values for which profile 
printouts are desired, m (ft) 

PRT 


Nt:, DEFAULT =0.95 

pr , t 

PRTAR 


N^ . array, used only if KODPRT = 3; 
Pr , t 

NUMBl values 

PTl 


p* , Pa (Ib/ft^) 

t ,oo 

R 


R* , gas constant (eq. (6)), m^/(s^-K) 

^ (ft^/ ( s 2-°R) ) ; DEFAULT = 1716 Et^/(s^-°R) 

SMXTR 


Critical-vorticity Reynolds number; 
DEFAULT = 1 X 10® 


30 


SST 


ORJQINAi: PmE FS 
OF POOR QUALITY 

S*-location of transition, s. . , m (ft); 
DEFAULT = 1 X IqS 


TLNGTH 

TTl 

VELEDG 

VISICI 

VIS1C2 

VIS2C1 

VIS2C2 

W 

WAVE 


i' transition extent (see fig. 2); 
DEFAULT =2.0 

T;,oc' K (°R) 

Value of F to be used in defining edge of 
boundary layer; DEFAULT = 0.995 

a* (see eq. (7a)), Pa-s (Ib-s/ft^) ; 

DEFAULT = 2.27 x Ib-s/ft^ 

a 2 (see eq. (7a)), K (^R) ; DEFAULT = 198.6 °R 

a^ (see eq. (7b)), Pa-s (Ib-s/ft^) 

a^ (see eq. (7b)) 

0 - neglect transverse curvature; 1 - include 
transverse curvature; DEFAULT = 0 


0 ^* 


^ is=0 
deg 


, shock-wave angle at s = 0 (see fig. 2), 


XEND 

XK 




k, constant in geometric progression (see eq. (42)) 


XMA 

M 

00 





XTl 


(see eq. 

(16a) ) ; 

DEFAULT = 

0.4 

XT 2 

”^2 

(see eq. 

(16b) ) ; 

DEFAULT = 

0.0168 

XT3 


(see eq. 

(17c)) ; 

DEFAULT = 

5.0 

XT4 

^4 

(see eq. 

(17c)) ; 

DEFAULT = 

0.78 

XT5 

k5 

(see eq. 

(20b) ) ; 

DEFAULT = 

0.108 

XT6 


damping 

function; 

DEFAULT = 

^ 26.0 


Input data for $NAM2 . 
Variable name 


Variable description 

Order of interpolation to be used for $NAM2 tables; 
DEFAULT = 1 


NUMBER 


Number of values read into $NAM2 tables 
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QW 


RMI 


RVWALD 


S 


SS 


TW 

Z 

Input data for $NAM3 . - 
Variable name 


Edge pressure-distribution array (NUMBER values ) , 
Pa (lb/ft2) 

Wall heat- transfer-distribution array (NUMBER 
values); KODWAL = 2 , W/m^ (Btu/ft^-s) 

Body radial-coordinate array r^ (NUMBER values), 
m (ft) 

Wall mass-flux array (NUMBER values) ; 

DEFAULT = 0, Pa-s/m (Ib-s/ft^) 

S-station array (NUMBER values) ; independent 
variable for tabular input, m (ft) 

Array of incremental values between adjacent solu- 
tion stations (s*, s*, s* ; step size 

1 2 lENDl 

can be arbitrary and not directly associated 

with the S-station array other than 

lENDl 

E SS (I) = the first two members of the 

max 

m=l 

array must be equal (SS(2) = SS(1)) 

Wall- temperature-distribution array (NUMBER 
values); KODWAL = 1, K (^R) 

Axial-coordinate array (NUMBER values) , m (ft) 


Variable description 


NUMBER Number of values read into $NAM3 tables 

RRS Array of radial coordinates of shock wave (NUMBER 

values) , m (ft) 

ZZS Array of axial coordinates of shock wave (NUMBER 

values) , m (ft) 

A unique relationship exists between the print control parameters (PROINC; 
PRNTINC; IPRO; IPRNT; PROVAL; PRNTVAL) and lENDl in $NAMl and the SS array in $NAM2. 
The potential user of program VGBLP should note that there are exactly lENDl values 
in the SS array and that these values specify the solution-station locations along the 
s-coordinate - Also, a solution station must be located at the s-coordinate locations 
designated as print (profile or wall) stations. A failure to understand this rela- 
tionship generally results in a computer run with no output. Consider the following 
input where the program user wishes to march the solution to = 1.0: 



SS = 10*. 001, 99*. 01 


ORIGfNAL PAGE !S 
OF POOR QUALITY 


lENDl = 75 

PRNTINC = 0.201, PROINC = 0.501 

IPRO = 1, IPRNT = 1 

PROVAL = 0.751, PRNTVAL = 0.751 

Two errors have been made in the preceding input that will result in the program 
stopping at s = 0.66 (instead of s = 1.0) without any output (wall print or profile 
print). The two errors are as follows; (1) lENDl is not equal to the number of 
values in the SS array; (2) the designated print locations do not agree with the 
solution stations designated by the SS array. An example of correct input is as 
follows : 


SS = 10*. 001, 99*. 01 
lEND =109 

PRNTINC = 0.2, PROINC =0.5, 

IPRO = 1, IPRNT = 1, 

PROVAL = 0.75, PRNTVAL =0.75 

The program would now have a normal STOP at s = 1.0 with wall prints at s = 0.2, 
0.4, 0.6, 0.75, 0.8, and 1.0 and profile prints at s = 0.5, 0.75, and 1.0. Finally, 
it should be noted that the SS array can be composed of completely arbitrary As 
values with the restriction SS(1) = SS(2), but to obtain output the user must specify 
print-control input corresponding to the location of the solution stations. 


Intermediate Data Storage 

The output (S, PE, RMI, TW, Z, DPEDS, RVWALD , DRDZ, QW) required at station m+1 
generated in subroutine TABLE is written on TAPE 4. Program VGBLP reads this output 
just prior to obtaining the boundary- layer solution at station m+1. For cases where 
variable entropy is included (lENTRO = 2), TAPE 4 is rewound at the end of the last 
computed station enable restart for the next variable-entropy iteration. 


Output Description 

Program VGBLP first prints namelist data for $NAM1. Next, subroutine TABLE 
prints $NAM2. If the case includes the effect of variable entropy (lENTRO = 2), sub- 
routine VARENT then prints $NAM3 input. It should be noted that for many cases the 
user can take advantage of many of the DEFAULT values for $NAM1 and $NAM2. 

Following the input data prints, the similar-solution profiles at the initial 
station (^ = 0) are printed as follows: 
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Variable name 



Free-stream and reference variables are then printed. 


AAl 

a*, m/s (ft/s) 
00 

PREF 

p*u* 2 , Pa (Ib/ft^) 

PTR 

p. /(p ) 

,00/ V'^00 00 / 

PTI 

p* , Pa (Ib/ft^) 

t , oo_ 

PI 

p*. Pa (Ib/ft^) 

00 

REY 


RREF 

p*, kg/m^ (Ib-s^/ft^) 

RTl 

Pt r kg/m^ (Ib-s^/ft^ 

t , 00 

Rl 

5 <g/ni^ (Ib-s^/ff^) 

TREF 

u*^/c*, K (^R) 
r / p 

TTl 

T* , K (^R) 
t,oo 
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T1 

T+ 

oo' 

K 

R) 

UREF 

* 

r 

m/s 

(ft/s) 

U1 


m/s 

(ft/s) 

VISREF 


Pa-s 

(Ib-s/ft^) 

XMA 

M 




oo 


The profile-print and wall-print stations are next printed in accordance with input 
specified in $NAM1 . 


Laminar-prof ile 
CROCCO 

ETA 


"'t - 

Tt,e - Tw 


FZ 

M/ME 

PT/PTR 

T/TE 

TT/TTE 

TZ 

U/UE 

VORTREY 

XLMll 

Y/YE 




m+1 f n 


M/M , Mach number ratio 
e 


p^/p^^^, total pressure ratio 


0 


V\,e 


90 \ 

9n/ 


m+1 , n 


^m+1 , n 


'^(hL\ 

V \9y / 


m+1 , n 


vorticity Reynolds number 


(PP) 


m+1 , n 


(PP), 


y/y. 


Additional values for transitional and turbulent profiles 


UDEF 

- ■ 


Uj 

UPLUS 

u 



VISEFF 

1 + ^ 


IJ 

YPLUS 

yu^ 
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Auxiliary-profile values (NAUXPRO = 1) 

1 - exp 


C 


DAMP 
EP 

EPl 

EP2 

FCl 

GRAD(U/UE) 
GRAD (T/TE) 

MIXDEL 


V 

Wall prin t 
BETA 

CFE 


K) 


m+1 , n 


^ /m+1 , n 




m+1 , n 



/p 

9u 


9y 

8f 


3n 


30 


8n 



(see eq . (16a) ) 


(see eq. (16b) ) 


/m+1 , n 


^ /m+1 , n 


V (see eq . (26) ) 

B, pressure-gradient parameter (see eqs . (30)) 

Tx, 


w 


1 2 
TT P U 
2 e e 

condi tion 


, skin-f riction coefficient based on edge 


CFW 


DLTAST 


DPEDS 


DSMXO 


DTEDS 


w 


, skin- friction coefficient based on wall 
density 

6*, displacement thickness, m (ft) 

dp. 


ds 

/3X. 


, pressure gradient 


max j 

Li 


ds ' 


temperature gradient 
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dUe 

velocity gradient 

value of DIF at convergence (see eq. (62)) 

A1 

0 


HD 

ITRO 

ME 

MUE 

NOITER 

NSTE 

NSTW 

NUE 

NUW 

OMEGA 

PE 

P20 

QSD 

RE 

REDELT 


RES 


% 

V heat- transfer coefficient, W/m -K 


T ~ T 
w aw 


(Btu/ft^-s-°R) 

number of iterations performed for variable entropy 
Mg f edge Mach number 
]Jg , Pa-s (Ib-s/ft^) 

number of iterations required for convergence 


^p(Pu)e 


•/ Stanton number based on edge condition 


c p u 
p^ w e 


/ Stanton number based on wall condition 


Nusselt number based on edge condition 
Nusselt number based on wall condition 


Pe' pressure. Pa (Ib/ft"^) 


^t,e 

2 

P u 

heat transfer, W/m^ (Btu/ft -s) 

P^, edge density, kg/m^ (Ib-s^/ft'^) 

p u 6* 
e e 

Po 


f Reynolds number based on local displacement 


thickness 


Pe^e^ 

— , local Reynolds number 

M tn 
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RETHET 


RFTRUE 

RMI 

ROUSE 

RSHK 

RVWAL 

RVWALD 

S 

SWANG 

TAUD 

TE 

THETA 

TRFCT 

TW/TTl 

UE 

UTAU 

VW 

XAL 

XT 

YE 

YMP 

Z 

ZSHK 


OF POOR QUAUit 


p u 0 

G 0 

, Reynolds number based on local momentum 

^e 

thickness 


T - T 

aw e ^ ^ 

, recovery factor 

T T 
t e 


body radius (see fig. 2 ), m (ft) 


X 

max 


local radius of shock wave, m (ft) 



(pv)^/(pu)^ 

(pv)^, dimensional mass flux at wall. Pa-s/m 
(Ib-s/f t^) 


s, boundary- layer coordinate (see fig. 2), m (ft) 
local shock-wave angle, deg 

2 

T^, wall shear stress. Pa (Ib/ft ) 

T^, edge temperature, K (^R) 

0, momentum thickness, m (ft) 

r, intermi ttency distribution (see eq. (38)) 


T 


w 


T 


t ,oo 


edge 



velocity, m/s 
m/s (ft/s) 


(ft/s) 


^m+l , 1 

(Y - 1)M^ 

6^, boundary- layer thickness, m (ft) 
n- value at y^ 

z, axial coordinate of body (see fig. 2) , m (ft) 
axial coordinate of shock wave, m (ft) 
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Flow charts and listings for program VGBLP and its subroutines are presented in 
appendix C. 


SAMPLE CASES 

A range of test cases is presented as guides for assisting users of program VGBLP 
in their own specific applications. Five major test cases are presented that include 
external and internal flows, flows with wall-mass transfer, flows where transverse- 
curvature effects are important, and flows where variable-entropy effects must be 
included. For each test case presented, the following information is given: (1) sche- 

matic of geometry; (2) boundary conditions; (3) all input data including variable 
dimension specification; (4) samples of output (see appendix D) ; and (5) plots of 
selected results. It is suggested that users compute two or more of the test cases 
prior to applying program VGBLP to their own particular problem. This approach is 
beneficial in that it (1) confirms that the software has been correctly implemented on 
the user’s computer system and (2) provide.s experience in using the program and speci- 
fying the correct input data; however, the user need not understand the algorithm in 
order to successfully apply program VGBLP. The first test case, flat-plate flows, is 
especially useful in developing experience with the grid specif ication and control. 

Test Case No. 1 

This case represents the simplest class of flow that is usually encountered. 

The flat-plate boundary layer need not be similar; for example, turbulent flow, arbi- 
trarily distributed wall-mass transfer, arbitrary heat transfer, or externally imposed 
pressure gradients result in nonsimilar boundary- layer development. 

For the present case, the test conditions of reference 22 are selected. A sche- 
matic of the model, including flow conditions and numerical results, is presented in 
figure 3. The input and sample output are presented in appendix D. Comparisons 



Figure 3.- Test case no. 1. 
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of numerical results with the experimental data of reference 22 are presented in 
reference 3. A grid refinement study showing the order of accuracy of the numerical 
approach is presented in reference 20. 


Test Case No. 2 

This test case is for the flow past a waisted-af terbody configuration (ref. 23). 
Transverse-curvature terms must be rncTuded; also, the flow is supersonic with an 
attached shock wave. A schematic of the model, flow conditions, and typical numerical 
results are presented in figure 4. The input and sample output are presented in 
appendix D. Comparisons of the numerical results with experimental data are pre- 
sented in reference 3. 



Figure 4.- Test case no. 2; skin- friction coefficient- 


The pressure distribution was taken directly from the experimental data (ref. 23). 
Results for two calculations are presented: (1) without transverse-curvature (TVC) 

terms; (2) with TVC terms. This particular configuration is an example of a body 
where the boundary- layer coordinate s cannot be expressed as an explicit function 
of the body-coordinate system R,Z, and as such must be obtained by numerical integra- 
tion. In the previous example for flat-plate flow, this presented no difficulty 
since S and 2 were congruent. It is suggested that the user of program VGBLP 
develop software to numerically generate the s-coordinate from a specified body- 
coordinate system and to interpolate edge-pressure and wall-boundary data, often 
specified as a function of the body-coordinate system, to the S, R-coordinate system. 
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tJ^spGzoidal rule is sufficiently accurate 
grate the following relationship: 


and can be easily implemented to inte- 



( 66 ) 


■rest Case No. 3 

Flows with wall-mass transfer are often encountered and can be efficiently 
solved by program VGBLP. The sample case selected is that of reference 24 for laminar 
boundary- layer flow. A schematic of the model, including flow conditions and numeri- 
pr'esented in figure 5. The required input and sample output are pre- 
sented in appendix D. It should be noted that program VGBLP can be applied to turbu- 
lent flow with wall-mass transfer if the user modifies the A+ definition in 
subroutine TURBLNT. (See ref. 2.) 


The ni^erical results for three wall-mass transfer boundary conditions are pre- 
sented in figure 5: (1) (pv)^<0 (suction) ; (2) (pv)^ = 0 (solid wall ); and 

U) (pv)„ > 0 (transpiration) . The input and sample output are presented in 
appendix D. Comparisons of the numerical results with experimental data are pre- 
sented in reference 3. It should be noted that program VGBLP is not limited to fixed 


values of (Pv)^; that is, (Pv)* = g(s) can be input in $NAM2 if required. 


Free-stream conciitions 


p. =4.14 MPa 
t ,«> 



(a) Skin-friction coefficient. 
Figure 5.- Test case no. 3. 
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(b) Velocity profiles; s = 0.100 ui. 


Figure 5.- Concluded, 


Test Case No. 4 

For hypersonic, blunt-body flows, the effect of variable entropy introduced by 
the bow shock wave can significantly affect the boundary-layer development. As an 
example, the flow over a 45° spherically blunted cone in helium flow is considered. 
(See ref. 25.) A schematic of the model, including flow conditions and numerical 
results, is presented in figure 6. Experimental pressure data, supplied by the 
authors of reference 25, were used as input. The shock-wave data were obtained from 
figure 7(d) of reference 25. The input and sample output are presented in appendix D 


Test Case No. 5 

Boundary- layer solutions are usually required for the design and analysis of 
nozzle flows. A typical wind-tunnel design case (see ref. 26) is presented in fig- 
ure 7. It should be noted that the solution for this case is initiated in the stagna- 
tion chamber and marched downstream to the nozzle exit. The input and sample output 
are presented in appendix D. Extensive comparisons of the numerical results with 
experimental data are presented in reference 26. 

Nozzle flows are typical of the more difficult applications of boundary-layer 
theory because of the large variation of pressure gradient dp^/ds as the solution 
proceeds from the settling chamber through the sonic throat and into the supersonic 
region of the nozzle. The thinning effect of the pressure gradient on the boundary- 
layer thickness in the throat region of the nozzle necessitates care in selecting the 
grid distribution in the normal-coordinate as well as the marching-coordinate 
direction . 
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A coiranon problem often encountered in large pressure-gradient flows is oscilla- 
tions in 6* caused by the specification of p^ = g(s) and/or step size in the 
s-coordinate - For many applications (e.g., rocket nozzle design), these oscillations 
may be acceptable, but for facility design, where the design goal is usually to 
achieve a shock-free flow that meets the design test-section flow conditions, caution 
and judgment must be exercised in specifying the inviscid pressure distribution and 
step-size distribution in the s-coordinate. For example, if a relatively course 
distribution of p^ = g(s) and a fine distribution of As were specified, the 

dpe 

resulting pressure-gradient distribution ~ • = 

functions for linear interpolation. Spline functions or higher-order interpolation 
could be used to obtain p^ and dp^/ds at the solution stations from the specified 
input values; however, care must be exercised since higher-order interpolation can 
introduce oscillations resulting in changes in the sign of the pressure gradient. 
Splines with tension represent the optimum technique for generating p^ and dp^/ds 
at the solution stations from the input data; however, experience has indicated that 
it is more efficient to input a sufficiently dense distribution of p^ = g(s) and 
use linear interpolation. It is suggested that the user of program VGBLP work several 
problems with large pressure variations in order to gain experience in specifying 
pressure input and step-size distributions. 


d[g(s) ] 


would be a series of step 


: 1 I 
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CONCLUDING REMARKS 


A computer program, VGBLP, has been presented for solving the compressible 
laminar, transitional, or turbulent boundary- layer equations for planar or axi sym- 
metric perfect-gas attached flows. A three-point implicit, variable-grid finite- 
difference procedure is used to solve the governing equations. The algorithm and 
software are modifications of the procedures presented in NASA TR R-368 and NASA 
TM X-2458, respectively. The modifications render the approach easier to implement 
while increasing the efficiency (computer resources) and accuracy as compared with 
the original approach presented in NASA TR R-368. 

Test cases have been presented and should serve as guides for potential users 
the software. These cases cover external and internal flows including flows with 
wall-mass transfer effects, transverse-curvature effects, and variable-entropy 
effects . 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
November 6, 1981 
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APPENDIX A 


DIFFERENCE RELATIONS 

Three-point implicit difference relations are used in references 3 and 4 to 
reduce the transformed momentum and energy equations (eqs. (28) and (29)) to finite- 
difference form. The differencing scheme proposed by Blottner (ref. 7) is used in 
program VGBLP. For completeness, both differencing techniques are presented. 

It is assumed that all data are known at the solution stations m-1 and m. 

(See fig. 1.) Then, it is possible to obtain the unknown quantities at the grid 
points for the m+1 station. In the subsequent development the notations G and H 
are utilized to represent any typical variable. 

Taylor-series expansions are first written about the unknown grid point (m+l,n) 
in the ^-direction as follows: 


m, n 


^m+l,n 


(G; 


V+l,n 


2 

AC2 


(Grr) 

m+1 , n 


3 

A^2 




(Ala) 


and 


1 „ = - (A^i + ACo) (Gr) ,, + 

m-l,n m+l,n ^1 ^2 t, m+l,n 


(A?^ + A^2>' 


(Grr) 


^ m+1 , n 


(A^, + A^^) ' 


'(Grrr) + - 

m+l,n 


(Alb) 


where subscript notation has been utilized to denote differentiation; for example, 
G^ = 3G/3^. 

Equations (Ala) and (Alb) can be solved to yield 




^1*^111+1,0 


X2Gm,n + Ag^ (A^ + A^s) 

2 Ag, 6 


^ggg 


and 


ACiA^2/ ^^2\ 

G = X.G - XcG . + 1 + TT— Grr + . . . 

m+l,n 4 m,n 5 m-l,n 2 \ ^^1/ 


(A2) 


(A3) 


Terms of the order of A^j ^^2' smaller, are neglected. This produces truncation 
errors of the order of AC]_ A ^2 instead of A^2 reference 5 where two-point 

difference relations are used. The X^^, X 2 r ••w coefficients appearing in 

equations (A2) and (A3) are defined as follows: 
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and 


+ 2 A? 


2 


ACi + A ?2 

A^l 


A?2 

^ AqcAq + AC 2 ) 

Aq 



(A4) 


(A5) 


(A6) 


(A7) 


(A8) 


Taylor series expansions are next written about the unknown grid point (m+l,n) 
in the n-direction as follows: 


APr 

®m+l,n+l ^m+l,n ^ ,n ^ 2 "'nn'm+l,n 


(G ) 

' irsTTi * 


An: 


■^^nnn^nHi,n + ■ • • 


(A9a) 


and 


3 


An, 


2 ^"^nn^m+l,n 


An 


n~I 

6 ^°nnn^m+i,n ^ 


(A9b) 


Equations (A9a) and (A9b) can be solved to yield 


3^2) ~ '^l‘^m+l,n+l " ^2*^111+1, n ^3‘^m+l,n-l 

/m+l,n 

(^’^n-l - An^) 

3 *^nnn ^ • 


(AlO) 
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and 

Arin 

'^4‘^m+l,n+l " ’^5‘^m+l,n " '^6^m+l,n-l 6 Snn + • • • 

(All) 


.^n/m+l , n 


The ^21 ■■■> Yg coefficients appearing in equations (AlO) and (All) are 

defined as follows: 


and 


Y-, = 


1 An (An + An -, ) 

n n n-1 


Yn = 


2 An„ An^_i 


■3 An„_i(An^ + An^_i) 


Ann_i 

""4 = An^(An^ + An~y 


^’3n-l - 


5 An^ An^_i 


Y. = 


An, 


Y^ = 


6 An^_;^(An^ + An^_i) 


(A12) 


(A13) 


(A14) 


(A15) 


(A16) 


(A17) 


tions 


For the case of equally spaced grid points 
{A4) to (A8) and (A12) to (Al7) reduce to 


in the and q-coordinates , 
the following relations: 


equa- 





(Al8a) 
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and 


= 


1 


An 


Y2 = 2Y^ 


Y. = Y, 
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(A18b) 


Y^ = 


4 2 An 


Y. = 0 


Y^ = Y, 


where A^ and An represent the spacing between the grid points in the and 
r)“ COO rdi nates , respectively. 

Equations (A2) , (A3), (AlO) , and (All) can then be written for constant grid- 

point spacing as follows: 


^3eLi,„ ■ 2 A? * 3 °KK 


(A19) 


m+1 , n in,n m-l,n ^ 


- 2G - G 


(A20) 




^m+l,n+l ^^m+l,n ^ ^m+l,n-l 


/m+l,n 


An^ 


12 nnnn 


+ . 


(A21) 


and 


9g\ ^m+l , n+1 ^in+l,n-l Ap^ _ 

^/mti,n ^ ~ ’ 6 nnn • 


(A22) 


Quantities of the 

linearized in order to 
this type are obtained 


form I^G -^j that appear in the governing equations must be 

obtain a system of linear-difference equations. Quantities 
from equations (A2) and (A3) . 


of 
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The procedure used to linearize nonlinear products such as 
same as that used by Flugge-Lotz and Blottner (ref, 5) and is as 



follows : 


is the 



(^] 

^'m, n ^m+1 , n 


\9n/m,n'^^/m,n 


+ 





(A23) 


where the terms 



and 


known station m. By equating 



are evaluated from equation (All) , but at the 
H in equation (A23) , the linearized form for 


quantities of the type 


3n/ 


is obtained; that is, 


. (^) 

V^^^m+l,n '^^^m,n 


2(-^'l - (- 


\9nym+l,n V 



+ 0(A^2)' 


where It:— I is obtained from equation (All), 

-^m+l , n 


(A24) 


The preceding relations for the difference quotients produce linear-difference 
equations when substituted into the governing differential equations (eqs . (45)), In 

references 3 and 4 it is noted that in practice the nonlinearities do not require 
iteration provided a sufficiently fine mesh-point distribution is chosen. However, if 
one wishes to increase the computational speed by reducing the number of grid points 
to a minimum, then iteration may become necessary . 


Blottner (ref. 7) proposed a variable grid scheme that in principle is only 
first-order accurate in the normal step size Aq^, but approaches second-order 
accuracy when the grid defined by equations (42) to (43) is used. In the variable 
grid scheme, the following difference relations are used: 


/^\ ^ ^m+l,n+I '^m+l,n-l 

\3n/ra+l,n '^^n-1 ^ 


2 




(A25) 



= Equation (AlO) 


m+1 , n 
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The nonlinear term is evaluated as follows: 




2 

1 1 
m+1 , n+— ' 

^‘^m+l,n+l ‘^m+l,n| 

L9n 

[ 3niJ 

m+l,n + An^.i 

V *"n / 


I / ‘^m+l,n ‘^m+l,n-l 

m+l,n— ^'^n-1 y 


i2|}‘^nnTi + ~ '^’^n-i^ 

where 

7 +7 

I - 1^+1/ n+1 ^m+l,n 

m+l,n+i 2 


(A26) 


(A27) 
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COEFFICIENTS FOR DIFFERENCE EQUATIONS 

Equations (45) are the difference equations used to represent the partial dif- 
ferential equations for the conservation of momentum and energy, respectively. These 
equations are repeated for convenience as follows: 

Momentum equation 


^^n^m+l,n-l ^ ^^n^m+l,n ^ ^^n^m+l,n+l ^ ^^n®m+l,n-l 


^^n®m+l,n ^ ^^n®m+l,n+l '^^n 


(Bl) 


Energy equation 


^^n^m+l,n-l ^ ^^n^m+l,n ^ ^^n^m+l,n+l ^^n^m+l,n-l 


^ ^^n®m+l,n ^^n®m+l,n+l ^^n 


(B2) 


These equations are obtained from equations (28) and (29) and the difference quotients 
presented in appendix A, The coefficients ^^n' forth, in equations (Bl) 

and (B2) are functions of known quantities evaluated at stations m and m-1. (See 
fig. 1.) Therefore, equations (Bl) and (B2) can be solved simultaneously. In refer- 
ences 3 and 4, equations (Bl) and (B2) were solved simultaneously without iteration. 
The reader interested in the coefficients for the noniterative simultaneous solution 
is referred to appendix B of reference 4. In the present approach, using the Blottner 
variable grid scheme (ref. 7), the coefficients are written as follows: 

Momentum equation 


^^m+1 , n^ q 

Al^ = - q 

n An^_l + An„ 




Y3. 


Bl^ = 


^ ^^m+1 , n^ g 




XI + 


(t^^Zs) + (f^^Ze) ,, 

m+l,n mtl ,n-H 


Y1 




Y3„ + 


(V 


mt 1 , n g 


An„_i + An^ 


ci„ = 


^^^V+l,n+l 


Yl. 


(B3a) 


(B3b) 


2 


(B3c) 
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Dl„ = 0 


El„ = -S 
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Pin = 0 




(X3) 




Energy equa tion 


2 (alt^^e] 


m+1 , n 


^2n = 


9f \ 

^E/iti+l , n 


A'^n-1 + A^n 


B2n = 0 


2 (aZt^^e) 


m+1 , n 


C2 = - 


^'^'m+l,n. 


Arin_i + APn 


(^m+l,n>q 

D2 = - ^ 

Ailn_i + APn 




Y3 


E2n = 


^ ^^m+1 , n^ 


AC, 


XI + 


(t^^Ze) . + (t23Zc) ,, , 

m+1 , n ^ m+1 , n-1 


Y1 


^^^^m+l,n m+1 , n- 1 


Y3_ 


^%+l,n^ 

^ ■ Arin_i + Ann 


F2„ = 


^ lu+l , n+1 

2 


Yl. 


^ ^^m+1 , ( 




[ix 2 ) 0 „_„ - (X3)0„.j_„] - („Zt2JF,^^nn|(S)" 


m+1 , n 


(B3d) 

(B3e) 

(B3f) 

(B3g) 

(B4a) 

(B4b) 

(B4c) 

(B4d) 

(B4e) 

(B4f) 

(B4g) 


These equations are solved using the coupled solution technique presented in refer- 
ences 3 and 4. A major difference between the present approach (program VGBLP) and 
that of reference 4 is that the system of equations can be iterated in the present 


53 



APPENDIX B 


approach as opposed to no iteration in reference 4. The iteration cycle is as 
follows: (1) the g subscripted quantities are first obtained from known values at 

stations m-1 and m using equation (A3) for the initial solution of the system at 
station m+1 (note that the continuity equation is integrated using the trapezoidal 
rule of integration) ; (2) for the remaining iterations the g subscripted quantities 

assume their value at the previous iteration at station m+1. The eddy viscosity is 
updated after each iteration. 


ORIGINAL PAGE IS 
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FLOW CHARTS AND PROGRAM LISTING 
Main Program VGBLP 


Program VGBLP controls the sequence of finite-difference solutions for the 
boundary- layer equations. It reads the input and, through its subroutines, sets up 
the computational grid, generates the initial solution, controls the parabolic march- 
ing procedure wherein the nonsimilar solutions are obtained, calculates all boundary- 
layer parameters, and prints the output at specified locations. The flow chart for 
the main program is as follows: 
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OF POOR QUALITY 

The program listing for program VGBLP is as follows: 


PkOGRAM VGBLP { in put# output^ TAPES-INPUT# TaPE6«UUTPUT> TAPE A) 
DIMENSION PTQJ>T(JK), STaB2(JK), MOMEIJK)# TTuTTUK)f CROCCO(JK), U 
lOUPL(JK)^ TCnRD(J<), U0EP(JK)» NuNDEL(JK), UEE(JH), PRDVAL(JN)# PR 
2 nTVAL(JM)» TAUP(JK)> RkS(JL)^ ZZS(JL)> DRSDZS(Jl), 0VT(JL>2), ANS( 
32)f PkTAR(JI) 

COMMUN /MESHl/ XN{JK)#DN(JK),Y1(JK)#Y2(JK)/Y3(JK),YA(JK)/Y5(JK)#Y6 
1( JK) 

COrtHON /SOL VI/ Fi(JK),F2(JK)/F3(JK),Ti(JK),T2(JK)#T3(JK),Vl(JK),V2 
1 (jK)»V3(JK,)»EP2(JK),EP3{JK)>FZ2(JK)»Pa3(JK)^T22(JK)>TZ3(JK),XL2(JK 
2)#Xl 3(JK),XLP2(JK)»XLP3(JK) > R A TOl ( JK ) » K AT02 ( J K ) > k AT 03 ( J K ) > E H2 ( J K ) , 
3EH3( JK)»DRAT01 { JK)>DkaTC 2 (jK ),DRAT03 ( JK)# VARA( JK)> VAKB( JK), VARC ( JK 
A)#VARO(J<)#VARE(JK)#Y(JKj#EPlUk)#EHl{JK)#XLl(JK)#XLPl(JK) 

COMMON /S0LV2/ Zi#Z2#Z3#Z<t#ZB#TE»FAA#FAa#FAC#FAD#BEXi#aEX2#BEX3,BE 
IXA 

EQUIVALENCE ( Z ZS # 0 V T < i# 1 ) ) # (ORSDZS>D VT (1,2 ) ) 

COMMON /TRBULNT/ S#KSTR#TLNGTH#Tki-ACT#0ISINC#XTl»XT2#XT6#XT3#XTA#X 
lTi.#PRrW,RE»UE#XN JE# J#RMI# EPS# jPOINT# IE#WWl#Wh2#WW3# WWA# WW5# NEOGE 
2# KOOViS# A# XBE » X# PR# KOu Pk f # Pk T# PRTAR# GLaR# NUMbl#XK 
COMMON /UNIT/ V I S 1C 1 # V I Si C2 # V I S2C 1 # VI S2C2# PTl# TTi # w A VE # R # PH 1 0#D S # S 
1ST#RT1# P1#TI#R1#U1# AAl# TREE# VISk£F#PESTAR#T£STAR#RESTAk#UESTAR#MU£ 
2STAR# YESTAR# THET A# TAUD# QSD# HD# U PLUS# D ISP# PE#2#Tw# QW#RVWALD# PROINC# 
3PRNTINC# ZS#RS 
INTEGER W 

REAL MOME#NJE#NJ»^#KrtD#KED# 1NTEGT#1NTEGL#N0NDEL#MUESTAK 
EXTERNAL INTEGT 

NAMELIST /NAMl/IGEQM# XEND#lE#XK#DETAi#XHA#PTl#TTl# IGAS# VIS1C1#VIS1 
1C2#VIS2C1#VIS2C2# G#R# PR#1 BODY# WAVE#PHI1# J#W# lENTRO# SST# SMXT R# KODV I 
2S#KTCiJD#TLNGTH# I Y I NT# KUDA MP # XT I# XT 2# XT3# X T A# XT 5# X T6# PR T # KOD PRT# NUM 
3Bi#6LAR#PRTAR# lENDl# PkClNC# PRNTINC#! PRO# 1PRNT#NAUXPR0#PR OVAL# PRnTV 
AAL# FT# KODE# KOD RAL# VELEDG# CDNV# N1TMAX#K0DUNIT# 1TMAX#C0NV E 

IN1TIALI2F DATA TO STANDARD INPUT 

DATA G/l.A/#R/17i6./#PR/. 72 / # PRT / , 95 / # w / 0/ # KQD E / 0/ # KOD W A L /I / # lENTR 
lU/l/#A/l./#K0DVIS/2/#SSl /i. EB/#FT/i. /#SMXTR/i.E8/»TLNGTH/2. /# XTl / 
2.A/#XT2/.0168/#XT3/5./#xTA/.76/#XT5/. 108/ # XT 6/ 26 . /# PRO INC/1. /#PRNT 
3INC/.l/#NAUXPR0/0/#NPUTYPt/l/# I PRC/ 0 / # i PR NT /O / # C ONV / . 000 1 /# NUMBl /O 
A/ 

DATA lHATS/0./#RS/0./#P20/0 ./#NiTMAX/l/#SI6N/l./#iT2lAX/i000001/#T 
1RFACT/0,/#KSTR/0/#ITMAX/3/» NOIT/0/ #0 PEDS/0, /# 2/0 , /# KOOPR T/1 /# KTC OD 
2/2/#KTCD/0/#IGAS/l/#K0DMMP/2/#VELEDG/.995/#IYlNT/i/#SMXN/0./#SMX0/ 
30./#SMXP/0./#?RTW/,95/#K0DUNIT/0/#ANRAD/1.7A5329E-2/#VISiCl/2.27E- 
Ad/#VISlC2/196,6/#XVALi/0./#C0NVt/.01/ 

Y(i)«0. 

NONOEL(1)-0, 

DO 1 I»1#JN 

1 Pk0VAL(I)«0. 

DC 2 I«1#JM 

2 PRNTVALI I)«0. 

DO 3 1«1#JL 
DkSD2S( I)*0. 
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3 UEE(I)«0. 

Du 4 I-1»JK 

EP2(X>*EP3(l)=XL2(l)=XL3(i)=KAT01tI)«RAT02(I)»RAT03(I)«Tl(I)«T2(I) 

l«T3(X)-Fl(I)»F2(I)«E3lI)*bPl(I)«1.0 

TAUP(X)«1. 

iTAd2(I)»F22(I)*rZ2{I)*XLP2(I)»XLP3{I)«FZ3(I)»TZ3lI)*DRAT01(I)«DRA 

XTD2 ( X }»DRATD3 ( I ) »0. 
i* CONTINUE 


READ NAMELIST INPUT 


READ (5,NAM1) 

IF (E0F(5)) 5^5 

5 STOP 3 

6 rR1TE16»82) IGEQM#XENU> IE#XR#DETA1#XMA#PTX/TT1/IGAS>\/IS1C1# V1S1C2» 

XVlS2CX»VIS2C2>G»R#PR>IBODt'>wAVE»PHlI»J^W#IENTRO>SST>SMXTR>KUDVIS#K 

2TC0U>UNGTH»IYINT»KUDAliP»XTi>XT2»XT3#XT<»#XT5#xT6#PRT»K0UPRT>NJMai# 

3GLAR/PRTAR 

RRlTE{6>83)IENDl»PKOiNC»PRNTXNC#iPRO/IPRNT»NAUXPRO#FT/KODE#KQDWAL# 
1VELEDG»CQNV^ NlTMAX>KQjUNlT» ITMAX»CONVE 
IF (KGDUNI T. NE . 1 ) GO TO b 
ATl-TTl 
AT2-TT1/2. 

IF ( IGAS.E0.2) GO TO 7 
AXX«vXSlCl+(ATl**l.ij)/{ATi + VISlC2) 
AX2»VIS1C1*(AT2*^1.5)/(AT2+VIS1C2) 

ACl»SuRT(5./9.)/^7.b7B2::3 

YISXCX»AC1*SQi<T(AXX*AX2’*‘(ATX + VISXC2)*(aT2+VIS1C2)/({ATX*AT2)^*1.5) 

1) 

VIS1C2*VIS1C2Y9. / 5. 

GO TO 8 

7 CONTINUE 

AX1»VIS2C1*(AT1**V1S2C2) 

AX2-VXS2C1*(AT2**VIS2C2) 

aC1-( (5./9, )■«<♦ VI S 2C2 ) /A7. 860 2 58 

VXS2C1»AC1 + SQRT(AX1*AX2/MATX*AT2) ♦♦ VIS2C2 ) ) 

8 CONTINUE 

P FII0«PHI I 

IF (KODUNIT.EQ.I ) CALL XnunXT ( PROV aL »PRNT VAL # JM# JN ) 

C 

C SET UP GRID normal TO wALL 

C 

Wl»XK 

W 2 ■ i ♦ m 1 

W3*I+wl+Wl^Wl 

RWl»W3 + W3 + (Wl*W2-l. ) + w2 

Ww2-wX*Vi2*W3*W3 

WW3»W3*W3 

Rv«A»l , +W1 

RV(5*ii«i + Wl*Wl*W2*43 

CALL MESH ( XK » X E ND^ IE » IGt OM » UE TAl ) 

RPR«1./PR 
DO 9 N-1# JK 
EH2(N)«RPR 
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EH3(N)«RPR 
EH1(N)-RPR 
9 CONTINUE 


C??!GIMAL PAGE f3 
OF POOR QUALITY 


PROGRAM CONSTANTS 

10 CONTINUE 

XMAC-1. + .5+ (G-1 . ) ‘^'XMA**2 
RTi-PTl/{R*TTl ) 

P1»PT1/(XHAC)++(G/(G-1. )) 

R1-RT1/(XMAC )*♦(!. /(G-1. ) ) 

TI-TTl/XMAC 
AAl-SORT (G^Pl/Rl ) 

Ul-XMA+AAl 

TREF»Ul+*2/( (G/(G-1. n+R) 

GPl«G-a. 

GMl-G-1 . 

GM10G*GM1/G 
XWAVE*ANRAD*WAVE 
IF (KODWAL.EO.i ) GO TO 12 
IF (XwAVE.EQ.O. ) GO TO 11 
XAF-(XMA*SIN( XWAVE ))**Z 

T20T1« (2.+G+XAF-GN1 )* (GHl+XAF+2. )/ {GP1*GP1*XAF ) 

XN2- (GPl+GPl+XHA ♦XMA+XAF-4. *(XAF-i . ) ♦(G+XAF+1, ) )/ ( ( 2 . *G *X AF -GMl ) ♦ { 
16M1+XAF+2. ) ) 

AWT»TI*T2aTl*(l.+SuRT(PK)+&Ml^XM2*0.5) 

60 TO 12 

11 AkT»TI+( 1.+SQRT( PR ) *G Ml +XMA*XHA*0. 5) 

12 CONTINUE 

IF (IGAS.60.2) GO TO 13 
VXSi-\/ISlCl*(TI^*i.5)/(Ti+VISiC2) 

VISRcF-VISlCl*(TREF*»1.5)/ ( TKEF+Vl S1C2) 

60 TO lA 

13 y/lSl»VlS2Cl*(TI**VlS2C2) 

VISREF»VIS2C1*(TREF**VIS2C2) 

lA CONTINUE 

REV-R1+U1+A/VI51 
REYREF-REY+VISl/VISREF 
XCaNE«ANRAD*PHIO 
EPS»1./S0RT (REYREF ) 

ABC-(XMA+SIN( XWAVE ) )**2 
P10»PT1/ (Rl + Ul + Ul ) 

IF ( Xki A VE.6T. .000 3001 . Ok . XW A V E . L T 0000001 ) P 10- ( 1 . / ( G»XMA ♦ X ,1A ) )♦ 

1( (XMAC+A8C* (G+1. ) ) / ( A6C*( G-i . )+2. ) )**(G/ (G-1. ) )♦( (G+1. ) /I2.*G*A8C- 
2( G-1. )))♦♦( 1 ./ (G-1. ) ) 

TlO-0. 5 + 1.0/ (XMA + XMa* (3-1.0) ) 

R10-6 + P10/ (T10*{G-i .3 ) ) 

TC-VIS1C2/ ( TI + XMAC ) 

RFL-SQRT(Pk) 

Rf-T-PK**0. 33333 
C 

C READ TABULAR 3ATA 

C 


IF (NOIT.GT.O) GD TO 20 
NA3C-NPUTYPE+1 
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GG TO (19»15>16)^ NABC 

15 CALL TABLE ( I E ND I ^ SO/ R 1 # Ui» A> TREF> KOD WAL# VI SR E K JDUNIT # AWT ) 
GO TO 17 

16 CALL TABLEl ( I E NDl/ SD> Kl » U1 / A/ TREF / K OOwAL / VI SR £ F/ KQDUN I T / A»< T ) 

17 GO TO (19/18)/ lENTRQ 

18 CALL VARENT ( R R S / Z Z S / D R i 0 Z S / NNN/ KOOU N IT / A ) 

19 CONTINUE 

SET UP PRINT control 


20 


SD-SD+A 

CALL SETUP (PROINC/ PROVAL/SO/ JN/IPRQ) 

Call setup (prntinc/PRntval/SD/JM/1prnt) 


.(RXT E(6/87) 

IF(KOOUNIT.EQ.l) WRIT£(6/8<t) 
IF(KOOUNIT.FO.l) WRiTt(6/85) 
1F(KD0UNIT.NE. 1) WKlTb(6/8A) 
IF(KODUNIT.NE .1 ) WRITt(6/65) 
CONTINUE 

KEAO (A) TW/RVWALD/GSO/PE/OS 


(PROVAL (I)^.30A6/ 1-1/ JN) 
(PRNTVAL(I)*.30A8/1-1/ JM) 
PROVAL 
PRNTVAL 


Si-DS 

UE-0. 

PEOPIO-PE/PIO 

IF { PtOPlO.GT. 1. ) GO TO 21 
GO TO 22 

21 CONTINUE 
ABTZ-ABSd.-PEOPlO) 

WRITE (6/81) ABYZ 
STOP 600 

22 CONTINUE 

IF (ibODY.NE.l) UE-SQkT(2.*T10»(l.-(PE0P10)**GM106) ) 

TE-TiO-.5+UE*UE 

RE»G*PE/{ (G-1. )YTE ) 

TR-TIO+TC/TE 

Tw-TW/TE 

SUT-VIS1C2/TREF 

IF (IGAS.EQ.l) XNUE-(TE**1.5)*((l.+SUT)/(TE+SuT)) 

IF (IGAS.EQ.2) XNUE-TE+YVIS2C2 
GO TO (23/25)/ I3GDY 


starting PRQCEOURE FOR FLOWS WITH STAGNATION POINT 


23 Sur-T10*TC 

IF (IGAS.EQ.l) VIS10»(T10Y*1. 5)^(1. +SuT)/(T10+SUT) 

IF (IGAS.FQ.2) VISiO»TiOY*ViS2C2 

P1PT2-P1/PT1 

If (XMA.LE.l. ) GO TO 24 

PlPT2»((((G + l.)*X'lA + *2)/2.)*Y(-G/(G-l.)))*((G + l.)/((2.*G*XMA**2)-( 
lG-1. )) )**{-!. /(G-l. ) ) 

24 CONTINUE 

JUEOS»SQRT(2.*((G-l.)/G)*TiO*(l.-PlPT2)) 

BELA-R10*VISI0*DUEDS 

RMI-SI 

0X10S»BELA+(SI-DS)YY(2*J+1) 


"I 1 1 
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DXDS-bELA+SI^*(2*J + U OF POOR QUALITY 

OXl-tBELA+SI+^CZtJ+ZII/IZ^J+a) 

X-OXl 

ARGH--778.26 + E PS*A*PR/ ( V i Sk E F *U 1 + U i + T 10 + SQR T ((J + 1>*R10*V1S10*0U60S 
1 ) ) 

XAL«0, 

XbE»l./ (FLOAT ( J ) 

GO TO 28 

starting procedure for flows without stagnation point 


25 continue 
IF (J.EO.O) GO TO 26 
RMI»SI*SIN(XCQNE ) 

BELEX«RE*UE*XNUE*(SiN(XCGNE )*+(2*J )) 

IF ( J.NE.O.AND.PHIO.EQ.O. ) 8 EL EX«R E + UE + XNUE 
DX10S-BELEX*((SI-DS)**(2*J)) 

GO TO 27 

26 CONTINUE 
RMI- 1 . 

BELEX-RE+UE+XNJE 
DXiUS'BELEX 

27 UXOS»3ELEX + 5I**(2*J ) 

OXl»(bELEX*(ST**(2*J+l)J)/(2*J+l) 

X-OXl 

QS-DISP»D. 

XAL»UE^*2/TE 
XBE«0. 

ARGM-0, 

DUEDS-0, 

28 CONTINUE 

GENERATE SELF SIMILAR SOLUTION 

CALL SIMILAR (I£» I G A S # X AL > X B E / PR # KOD WAL^ T W# D T D Z W> XK # TR# \/ 1 S2C 2» F 3 > F 
12»T3»T2»V3»V2^FZ3»TZ3>XL3#EP3»ARGM»OSO) 

IF (KODWAL.EQ.l. AND.lbODY.EQ.i) QSO« TZ3 ( 1 ) +XL 3 ( i ) / ARGM 
MUESTAR«XNUE*VISREF 
UESTAR-UE+Ul 
TESTAR»TE*TREF 
PESTAR»PE*Rl+Ul*Ui 
RESTAR-RE+Rl 
X*0. 

DX2-DX1 
C 

C SET UP INITIAL PROFILES 

C 

00 29 N-1»IE 
\/2(N)«V3(N) 

V1(N)«V3(N) 

T2(N)-T3(N) 

T1 (N)«T3(N) 

F2(N)«F3(N) 

F1(N)-F3{N) 
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XL2(N)»XL3(N) 

XLi (N)»XL3 (N ) 

XLP2(N)»XLP3(N) 

X L PU N) *XLP3 ( N ) 

TZ2(N)»TZ3(N) 

fZ2CN)«FZ3(N) 

29 CONTINUE 

IF (KOOUNIT.EQ.l ) CALL OuTuNlT ( PROV AL# PRNT V AL# J J N ) 

IF (KOOUNIT.EQ.l) call wALLUUT ( PkuV AL^ PR NTVAL » JMf J N ) 

WRITE initial STATION PARAMETERS 

WRITE (6#72) MUESTAR»UESTAK»TESTAR#PESTAR^QSD#XAL»X8E/RESTAR 
If (KODUNIT .NE . 1 ) WRITE ib,o9) 

IF (KOOUNIT.EQ.l ) WRITE (5»70) 

PTREF-PTl/ (R 1*'J1 + U1 ) 

WRITE (6»71) PTl*ni>RTl>Pl>ri#Kl^Ul»AAl#XMA>KEY#kl+Ui+Ul#TREF#Rl# 
1U1»\/ISPEF»PTREF»P10>T10/k1C 
IF (KOOUNIT.EQ.l) CALL INUNIT ( PRO V A L # PRNT V AL # J M# JN ) 

ZS-0. 

IF( (1ENTRD.EQ.2) .AND. (NOIT.EU.O)) WRITE(6#68) 

3EGIN NARCHING PROCEDURE ALONG SURFACE 
SMI— OS 

S ■ 0 . 

0XDS*OX10S 

IENDPi-IENDl+1 

00 68 M-2» lENOPl 

SM2*SM1 

Shl»S 

READ (A) S» PE^RMI .TW>Z^DPEDS»RVwAlO»DRDZ#QW 
PHI-aTAN( DROZ ) 

C0STH*C0S (PHI ) 

PP*DPEDS 

jE«SQRT(2.0^T10*(1.0-(PE/P10)**( (G-1 .0)/G) )) 

IF (NUIT.GT.O) J£*UEE(M) 

TE-T10-0.5+UE+UE 

XAL-Ue+UE/TE 

RE»G*PE/((G-1.0)^TE) 

IF (IGAS.EQ.l) XNUE» ( TE++1. 5 )♦ (1.0+T 10+TC )/ (TE+TIO+TC ) 

IF (IGAS.EQ.2) XNUE«TE**V1S2C2 

0X20S-DX1DS 

OXIDS-OXOS 

0X0S-RE*UE*XNUE*(RMi'i‘*(2*J)) 

DX2«(DXDS+DXlD3)*.5*Di 
IF (I80DY.EQ.1) OX2»OX2*.5 
IF (M.EQ.2) 0X1-DX2 

IF (N.EQ.2) GO TO 31 

CHECK FOR CHANGE IN STEP INCREMENT 

CKK-IS-SMD/ISHl-SMZ) 

IF (CKK.LT. .999999999 ) Gu TO 30 
IF (CKK.GT.i.OOOOOUi) GO TO 3o 
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0X2»(S-SM1) + {DXOS+4.*axlOS-t-DX20S)/3.-DXl 
GQ TO 31 

30 0X2»(S-S«l)*{DXlt)S+DxDi )/2. 

31 X-X+DX2 
OY»2.^X 
OZ»SURT(Oy) 

DUEDS*-PP/(»E*UE*OADb) 

0TEUS»-UE*DUEDS 

XAL-UE*UE/TE 

XbE-OY+DUEDS/UE 

Tk-TiO*TC/TE 

SET UP OIFFERENCE-OUOTI ENT CUEEFICIENTS FOR X-CJOROINATE 

Zl»2.^((DXl+2.*DX2)/(jXl+UX2)) 

Z2*2.*OXl + DX2 )/DXl 
Z3-2.+((DX2*0X2)/(DXl*{UXl+DX2))) 

ZA«{ DX1+DX2) /DXl 
Z5-0X2/DX1 

FAA*OZ/(RE*UE*(RMI**J) ) 

FA3»2.*EPS‘»‘W*FAAkCaSrH/('<MI**J) 

FAC-2.*EPS + 0Z=»'C0STH/(KE»UE + RM1+RM1) 

FAD»RMI/(EPS+C0STH) 
ttEXl»(Rm* + J)/(EPS*CaSTH) 

IF (J.EQ.O) SIGN«ABS(3IGN) 

BEX2-SIGN** J 

BEX3«(2. **JM'EPS*COSTH«SQkT (2. +X) /{R E+UE* (RMI ♦♦ (2. + J) ) ) 

BEX4-1. / ( J*1 ) 

DFOZWl. 

NlT-0 

32 CD.>JT1NUE 
NIT-NlT+1 

IF (KGDWAL.NE .1 ) AR GH = - QW *7 76 . 2 &♦ ( E P S+A / ( V I SF E F»U 1**2 ) ) ♦ ( PR *OZ/ ( R E 
1*UE»TE*XNUE*RMI*^J ) ) 

CALL SOLVE ( KOD WAL # AkGH^ T R, V IS2C2 , XAL # XBE ^ X# DX2# T IGAS> IE# J*NIT>R 
lVkALD#Rl»Jl#RRl»RW2>WW3>RwA^W<*i>,XNuE^XK#EPS>M,IB00Y>FT) 

C0»0. 

DiSINC-0. 

0ISP«0. 

TPK-0. 

TFIE TA*0 . 

KON-IE+2 
OU 33 N»2>KDN 

TPK»TPK+,5*{T3(N-1 )+T3(N) )*DN(N-1) 

RAT03(N)«SIGN*SQRTa.+FA3*TPK) 

0RAT03(N)-FA3*T3 (N) 

Y(N)»BEX1>M-1.+3EX2*(1. + BEX31‘TPK)**BEX^) 

DISP«DISP-»-.5*{(T3(N-l )-f 3 (N-1) ) + (T3 (n)-F 3(N) ) )^DN(N-i) 
C-1./(RATD3{N)*+J ) 

C»C + (F3{N)t(l.-F3(N) } ) 

THETA«THETA-*-.5 + (C0 + C )*ON(N-l ) 

CO-C 

D1SINC-DIS1NC+.5+1 ( T3(N-i )i( l.-F3(N-l))/RAT03(N-l))+(T3(N)* (l.-F3( 
IN) )/RaT 03(N) ) )+DN(N-l ) 
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33 continue 

0ISINC«0IS1NC*EPS*EAA*& 

THETA»THETA*EPS*F4a+A 

01SP«0ISP*EPS*FAAtA 

C 

C 

OETEkMINE THE LOCATION OF THE BOuNOaRY-LAYER EDGE 
C SEE INPUT DESCRIPTION FOR DEFINITION FOR VElEDG 

C 

NCOuNT-0 
DO 3A N»l>IE 
NCOUNT-NCOUNT + 1 
IF (F3(N) .GE. VELEDG) GO TO 3D 
3A CONTINUE 

35 NEDGE-NCOUNT 

Y£OGE-Y(SFDGE-l ) + ( VELEDG-F3 (NEDGE-1) )*(Y(NEOGE )-Y(NEOGE-D) / (F3(NE 
lDGE)-F3(NEDGE-i) ) 

CHECK THE CONVERGENCE 

DIF-At}S(l.-FZ3(l)/DFDZw) 

OFOZw-FZ3(l) 

DIFi«DIF 

IF (NIT.GT.NITNAX ) 0IF»0. 

CALCULATE WALL AND INITIAL VALUES REQUIRED FOR BASIC BOUNDARY 

layer parameters 

liOME(l)»F3 ( 1)/S0RT(T3(D) 

XMF1»2.+X4L 
XhF2»T3(l)*TE/T10 
XMF3-I.-XMF2 

Tr0rT(l)»(2.*T3(i)+XAL*Fi(l)*F3(l) )/ XMFl 
CkQCCO(l)»(TTQTT( 1 )-XMF2)/XMF3 
U0UPL(1)«0. 

TC0RD(1)»0. 

UD£F(1)*0. 

NCNOtL( 1 )»0. 

XMFA»EPS*XNUE*UF<‘UE*(RMI**J)*XL3(l)*FZ3(i)/QZ 

CALCULATE BASIC BOUNDARY LAYER PARAMETERS 

KON»IE+2 
DO 37 N»2>KON 

IF (TkFACT.EQ.O.O) GO TO 36 
UPlJS»SQRT(XMFA*T3( 1) ) 

UOUPL (N)»'JE + F3(N)/JPLJS 

TCORD(N)« ( Y{N) ^-UPLUS+RE ) / (EPS*XNUE*( ( T3 ( N ) + XL 3 ( N ) ) ♦♦2 ) ) 
UDEF(n)»UE*{1.-F3(N))/uPlUS 

36 CONTINUE 

rTOTT(N)»l2.RT3(N)+XAL*F3lNJ*F3(N))/XMFl 
CkOCCO(N) » (TTOTT (N)-XMF2)/XMF3 
MOME{N)*F3(N)/SQRT(T3(N)> 

37 CONTINUE 
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calculate VORTICITr RLYNGLOS NUMBER 
AND determine stability index 


IP ( TkFACT. GT. 0. 0 .AND . TKFaC T.LT. 0.9999) GO TO AO 
00 38 I«1>NEDGS 

TE«iia»RE+«E*UEtUE + Y{I) + Y(I)*(KMI*RAT03(I))* + J/(S&RT(2.*X)*XNUE*EPS 
1 ) 

38 5TA32(I )»Tf RM3+FZ3C I)/ (XL3( I )vT3(I )*+3) 

ST2MAX-STAB2 (1 ) 

OD 39 IX-2>IE 

IF (ST2MAX.GT.STA32UX) ) GO TO 39 
ST2MAX»STA0?(IX) 

39 CONTINUE 
ShXP»ST2MAX 

AO CONTINUE 

DSrtXO»{ SMXP-S.MXN) / ( S-SMl) 

SMXN-SMXO 

SMXO-SMXP 

IF (KTCOO. EQ. 2. GR. kTCD. Ew. 1 ) Gu TO A1 
KTCO-1 

RRRR»RE*UE*R1’«'U1 / ( XNUE + Vi GREF ) 

rLNGTH*l. + ( 'i. + RRRRA* (-1 )♦ (RRkk*A*SST )<‘*,8 )/SST 
A1 continue 

IF (SMXTR.LF.ST2MAX) CALL TURBLNT ( T 3> XL3# F Z3* R aT 03, Y> E P3> F 3> Eri 3# V 
lARA, VARB, VARC» VARD,VARE,.<Ul^AMP,IYINT,YEDGE,TAOP ) 

IF (SST.LE.S+A) call TjRBLNT ( T 3, X L 3 , F Z3» R A TO 3 , Y, t P 3, F 3 , EH3 , V AR A , V 
lARB, YARC»\/ARDfy/ARE,KaL)ANP,I YInT, Y£DGE,TAjP) 

IF (OIF.GT.CONV) GO Tu 32 

OUTPUT QUANTITIES 

DO A2 N»2,K0N 
NOnOEL(N)-Y(N)/YEDGE 
A2 continue 

IF (J.NE.l) GG TG AA 

TIBUR0N«2.*DISP/RMI 

IF( riBURON.LT.-l . )G0 TO A3 

DISP-RMI*{-1,+S0RT{1.+TI6UR0N)) 

GO TO AA 
A3 nRITE(6,93) 

AA CONTINUE 

THAi)lS»DISP/THETA 

REOELT-(RE+UE*OI5P/(XNUE^A))+REYkEF 
RETHET*REDELT/THADiS 
RES»(RE*UE*S/XNUS)*kEYREF 
XMAE-OE/SORT (TE*(G-1 . ) ) 

TW0TTi*T3(l ) +TE+TREF/ TTl 
RFCT0R*(T3(l)-i, ) / ( ( T T 1 / ( TE 1' TR £F ) ) -1 . ) 

TAUO»vISREF*Ul+ ( XL3 (1 ) »R E *X NUE ♦UE*UE ♦ (kMI J ) +FZ 3 ( 1 ) / OZ ) /(E PS*A ) 
CFE»TaUD/( . S*RIAJl*UH'Rt’PUE*Ut) 

CFWCFE + TW/TE 
IF (KuOWAL.NE .1 ) GO T J A5 

US»-XL3 ( l)*RE*JE*Tt + XNUE^(RMI-l‘*J) + TZ3{i)/(FR*0Z*EPG) 
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QS0«CiS*VI?REFtul*Ul/(77e.26*A) 

GU TD Afc 
A5 QSJ»QW 
f^b CDNTiNUF 

1 1- (FZ3(1) .GE .0. ) GO T O 92 
WRUE(6»91)S»T4UDfQSD 
STOP hOO 
92 CCMTiNUE 

Twa-TREF^TW 

TAwO-IRFL + (RFT-RFl ) MRFAC T) + (TTl-TE*TREF)-*-Tt*TREF 
Ht*OSD/{ TikD~TAWD) 

CMc«773,26*HD-*‘(G-l.)/{G*kt‘ftl*Ul*RE*UE) 

CHrt*CHE+TW/TF 
HSD»HD*A*S 

K£D»K*VISREF<‘XNJE*G/(7 78.26vP(;;*(G-1.)) 

KWD»KED*XL3(1)+T3{1) 

NUc-HSD/KED 
NU»<»HSD/KWD 

TOTAL PRESSURE RATIO# ( P T2 ) B L / ( PT2 ) R E F 

DO 49 I«1»IE 
ZEd»(HOME(I )^XMAE 
IF (SQRT(ZEP)-I.) 47#47#4B 

47 PTZ-PE^C (l.+(G-l. )*ZEb/2. )♦♦(&/( G-1. ) )) 

GO TO 49 

48 BBA»lZEB*(G+l.)/2. 

Bfa3»((G + l.)/((2.*G*ZEB)-(G-l. )))*♦(!. /I G-1.)) 

PT2»PE*BBA*6BB 

49 PTJPT(I)«PT2/PTR£F 

VARIABLE ENTROPY CALCULATIONS 

GO TO (5l»50), lENTRD 

50 lNT6Gt.»INTEGT(F3#NtOGE#DN) 

X X J* J 

RS»((J + 1)’*‘EPS*0Z»INTE&L)Y’*'((2.-XXJ)/2.) 

IPT»-i 

CALL lUNi { JL » NNN» RRS# 2# OVT # 2# kS# ANS# IPT# lERR ) 

ZS»ANS(1) 

TANThTS»ANS(2) 
that:>»atanitanthts ) 

SWANG-THATS/. 0174533 
PTTWO*(G-H. )X‘XNA + *2*bIN(THATS)X‘*2 

PTTkO»PTTwO/(IG-l.)*(XKA)*’F2*SlN(THATS)**2 + 2.) 
PTTkU»PTTWD’F*(G/(G-l. ))*PT1 
PTTkQ»PTTwO*{(3+l.)/ ( 2 . *G *X M A**2*S IN ( THAT S ) **2- ( G-1 . ) ) )♦♦(! ./ (G-1. 
1 ) ) 

P20«PTTk0/ (PI* 

UbE(N)-SQRT(2.*T10»(l.-(PE/P20)**((G-l.)/G))) 

XVAL2»ABS( (UEE(N)-OE)/UE) 

IF(XVAL2.GE.XVAL1 ) XNX-XVAL2 

XV'ALi»XVAL2 

NN-IE+2 
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51 COMTINUE 

IF (lENTRO.EQ.l) PZ0»P10 

PRINT PROFILES AND hALL VALUED 
KPRNT-0 

IF (ABS(A-1. ) .LE. 0.0000001) GO TO 52 

Z»Z*A 

RM»RMI+A 

ZS*Zb+A 

RS»RS*A 

52 CONTINUE 

DO 53 NUMBR«1#JN 

IF (S.GT. PRO VAL (MUMBR )-. 000001 .AND .S .LT .PRO VAL (NUMB R) + . 000001 ) GO 
ITO 5A 

53 CONTINUE 
GO TO 63 

54 CONTINUE 

IF (KODUNIT.EO.l ) S»S*.3046 
WRITE (6#73) 5 
IF (KODUNIT.EQ.l ) S»S/.304b 
I£OGEX»NEDGE+lO 
IF (lEDGEX.GT. IE) iEDGEX*lE 
IF (NAUXPRO.NE . 1 ) GO TO 56 
WRITE (6»74) 

DO 55 I»1,IEDGEX 

55 WRITE (6,75) ( X N ( I ) , T A J P ( i ) , V 3 ( I ) , F Z 3 ( I ) , TZ 3 ( I ) , V AR A ( I ) , \/ AR B ( I ) , V A 
lRL(I),VARD(I),EP3(i),VARt(l)) 

56 CONTINUE 

IF (^uDE.EQ.l) GO TO 5b 
IF (TRFACT.GT .0.7994) GU TO 61 

IF (TRFACT.GT .0.0. AND. TkFACT.lt. 0.9999) GO TO 58 
WRITE (6,76) 

DO 57 I»1,IEDGEX 

57 WRITE (6,77) X N ( I ) , NONDE L ( 1 ) , F 3 ( I ) , T 3 ( 1 ) , I TOT T ( 1 ) , C R JCC G ( I ) , P TO P T ( 
II ), NOME (I),FZ3(I),TZ3(I),STaB2(I),XL3(1) 

GO TO 63 

58 WRITE (6,76) 

00 59 I-1,IEDGEX 

59 WRITE (6,77) X N ( I ) , NGND £ L ( I ) , F 3 ( I ) , T 3 ( i ) , TTGT T ( I ) , C ROC C 0 ( I ) , P TO P T ( 
1I),M0ME(I),FZ3(I ),TZ3(I),6TA82(I),XL3(1 ) 

WRITE (6,76) 

00 60 I»1,IFDGEX 

60 WRITE (6,77) X N ( I ) , NGND £ L ( I ) , F 3 ( i ) , T 3 ( I ) , TT QT T ( 1 ) , C R ICC 0 ( I ) , PTO? T ( 
II),iiUME(I),TCORD(I ),g0UPL ( I ) , UDEF ( I ) , EP 3 ( I ) 

GO TO 63 

61 WRITE (6,78) 

DO o2 I«1,IE0GEX 

62 WRITE (6,77) X N ( I ) , NONDE L ( I ) , h 3 ( I ) , T 3 (1 ) , 1 1 QTT ( I ) , C ROCC 0 ( 1 ) , PTOP T ( 
II ),NOME ( I),TCQRD( I ),DOUPL ( I ) ,UD£F( I) , EP3( I ) 
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63 CONTINUF 

DO 6A NUM3ER«1»JH 

IF (S.GT.PRKTVAL ( NUMBER )-. 000001. AND .S.LT. PkNTVALCNUrtB E k)+. 000001) 
1 GO TO 65 
6^ CONTINUE 
GO TO 66 

65 CONTINJE 
PESTAR-PE+Rl+Ul^ J1 
TESTAR-TE*TREF 
RESTAR»RE*R1 
UESTAk-UE^Ul 
MUeSTAR-XNUE+VISREF 
rb3TAk»YEDGE+EPS*A 
RVWAL»RVWALO/(RESTAR*UESTAk) 

GANODG-RVWALD 

I E (KDUUNIT. EQ. 1 ) G AND uG» GA NDOG/ . 00636 5 B80A 
IF IKDDUNIT.EO.l ) CALL WALLOUT ( PROVAL# PRNTVAL# JM» JN ) 

WRITE (6»79) S>REThET#PP»CFw#ZS#YE 3TAR#X#RES»DTEDS#0SD#RS#UPLUS»RM 
iI#PE3TAR»DUEOS>HD»NQir»lRFACr>Z»TESrAR»DI3F#CHE#TwOTTl#JPOlNT^XBE# 
2RESTAR#THETA»CHW»RECTDK>P20 

WRITE (6»80) XAL^ UESTaR,THADaS,NUE,ST2P.AX#EP3#GAND0G#XMAE»TAUD#NUW 
1^ D3i1X0»REDE LT.MUESTAR>CFE» SwANG» V3(l )»RVwAL»NlT»DIFl 
KPRNT-1 

66 CONTINUF. 

AiNY*l./A 

IF (KODUNIT. EO.l .A nD.KPRnT.EU.1) A IN V»A I N V + 3 , 2 808 39 895 

3-i*AINV 

Z«Z*AINV 

RMl»R«I*AIN\/ 

Zi»ZS*AINV 

KB-RS+AINV 

UPDATE VARIABLES FUR MARCHING PROCEDURE 

DU 67 N«1 > NN 
F1(N)»F2(N) 

Ti(N)«T2{N) 

VI (N)*V2 (N) 

RaTU1{N)-RAT02 (N) 

JkATUl(N)»DRATQ2 (N) 

XL1(N)-XL2(N) 

XLP1(N)«XLP?(N) 

F2(N)»F31N) 

T2 (N)»T3(N) 

V2{N)*V3(N) 

RaTj2IN)«RATG3(N) 

DRATG2(N)»DRAT03(N) 

XL2(N)-XL3(N) 

XLP2(N)«XLP3(N) 

TZ2(N)»TZ3(N) 

EP1(N)«EP2 (N) 

EP2(N)*EP3(N) 

EPSM-EPSIN) 

IF (N.GT.l.AND.N.LT.iE ) E P3 M- ( EP3 ( N-1 ) + EP 3 (N ) +E P3 (N + 1 ) ) / 3 . 
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TAUP{N)«l./(XL3(l)*l-Z3{l))*XL3{N)*ABS(HZ3{N))*EP3i1 
£H1(N)«EH2{N ) 

EH2(N)*EH3 (N ) 
pZ2(N)«FZ3(N) 


67 CONTiNJE 

68 DX1»DX2 

IE (lENTRO.EQ.l ) STOP 100 
lE(NOlT.EO.ITMAX) i»RlTE(b/b9) 
IF (NOIT.FQ.I FMAX) STUP 200 
IFCXMX.LE.COMVE ) wRlTE(b>90) 
IF (XMX.LE.CONVE ) STOP 300 
NOir»NOIT+l 

WPITE(6>86) ITiAX>N01T 
REWIND ‘f 
GO TO 10 


69 FORMAT ( IHl, 20X» 5AH++++01 MENSiONAL OUTPUT QUANTITIES ARE IN U,S. S 
ITANDARD UNITS****W) 

70 FORMAT ( IHl , 20X, 5 5rl»t + + 0 i ME N S I UNAL OUTPUT QUANTITIES ARE IN S.l. U 
1NITS^***> / ) 

71 FORMAT (2X»30HFREE STREAM V ALUt S-D IM ENSI JNAL 5X> 5 HPT 1 «#E1<».6*2X 

i>5HTTl E1^.6f 2X» 5HRT1 » » E i ^ . 6> 2X# 5 H PI El<*. 6> 2X» 7H Tl-iE14. 
26/2X/5H R1 EU.6W> 5XoH Ui E 1 4 . 6> 2X» 5HA A I E 1 4. 6» 2X> 5HXMA 

3E14.6^2X, 5HREY - » E i 4 . 6 W / > 2 X , 2 8 HR E F E R EnC E V AL UE S-01 MENS I ONA L W » 5 X , 
55HPREF«,E14.6.2X,5HrREF*»E14.6^2X#5HRREF«»E14.6,2X/ 5HUR EF-# E 1 4 . 6, 2 
6X#7HVISREF»^ E14.6#2XoHPrK *>E14.6# / /2X# 2 5HPARAME TERb-NONOI MENS I ON 
7ALW>!?X,5HP10 =• f E 1 4 , 6» 2 X» bH T 10 » > E 1 4 . 6» 2X # 5riR 10 »»EI4.6W) 

72 FORMAT (2X»5HMUE » » £ 1 2 . o» 2x » 6HUE -» E 12 ,6» 2X » 5HTE Ei 2 .6, 2X>5HPE 

1 »#E12.6#2X»5HQSJ « # E 1 2 . 6/ / # 2 X> 5HX A L »# E 12 . 6# 2 XoHXBE */E12.6#2Xf 
25HRE »>E12.6) 

73 FORMAT ( / / 1 2 HS « F 14. 4# 9H PROFILE/) 

74 FORMAT (130H ETA TAUP V GRAO(U/UE) GRAD 

KT/TE) FCl DAMP EPI EPO EP 

2 MIXDEL/) 

75 FORMAT (11E12.3) 

76 FORMAT (//4X , 3HE T A > 8X > 4 Hr / Y E , 7X# 4HU/ U£» 7X # 4HT / T E> 6X>6HTT/TTEf 5X# 6H 
lCR0CC0»5X»6HPT/PTR^6X»4hM/Nfc»8X#2HFZ»9X»2HTZ>6X#7HVGRTREY#6X#5HXLM 
211/) 

77 FORMAT ( 12E11 .3) 

78 FORMAT ( / / 4X / 3HF T A# 8 X , 4HY / Y E> 7X, 4HU/ U£^ 7X, 4HT / TE# 6X , 6HT T / TT E > 5X , 6H 
1CR0CC0> 5X/6HPT / PFR# 6X, 4HM/ME WX, 5H YPLUS# 6x# 5HUPLUS# 6X# 4HuDE F > 6X, 6H 
2y/iScFF/) 

79 FORMAT (/2X#7HS E 1 2 . 5 # 2X# 7HRETHE T«# E 12 . 5# 2X, 7H0PE DS «^E12.5> 

12X,7HCFW *»E12,5,2X>7HzSHK El 2 . 5 #2X, 7HYE ■ # E12 . 5 W2 X, 7H X I 

2 -,tl2.5,2X, 7HRES - / E 12 . 5, 2X, 7H0 T EDS E12 . 5, 2X» 7HQSD ■»E12. 

35#2X#7HRSHK » » E 1 2 . 5/ 2 X , 7 HU T AU El 2 . 5» / 2X, 7HRMI E 12 . 5» 2X> 7HP 

4E »>E12.5#2X, 7H0UEDS = ^ E 12 . 5# 2X, 7HHD E 12 . 5# 2X> 7H I TR 0 »#I1 

52#2Xj»7HTRFCT « » E I 2 . 5 W2 X» 7H Z E 1 2. 5# 2X> 7HTE E 1 2 . 5 / 2 X# 7HD 

6LTaST-»E12.5>2X>7HNSTE * > E 1 2 . 5# 2X# 7HT w/ T T 1« > E 12 . 5# 2X# 7HY MP ->I1 
72>/2X,7HBETA E12 . 5^ 2X, 7HRE « # E 12. 5» 2X , 7HTHE TA ■ > E 1 2 . 5 » 2X, 7HN 

8STw »»E12.5^ 2X» 7HRFTRUE«» El2,5i2Xi7HP20 »>E12.5) 
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bO FORMAT (2X,7 HXAl « » E 1 2 . 5 » 2X WHUF » > E 12 . 1> # 2X, 7HF QRM -#E12.5#2 

1X,7HNUE «#E12. 5»2a# 7HFUUSE « # E 12 . 5 » 2X, 7HLME GA » , E 12 . 5 W 2X # 7HR i/WA 
2lD-,fcl2. 5» 2X, 7HME » > El 2 . 5» 2 X» 7HTA uO « » E 12 . i?# 2X> 7HNUW «#E12.5 
3»2X»7H0SMX0 «^E12.:j/2X,7HRE0ELT = #tl2.5,2X,7HMUE -,E12.5,2X,7HCFE 

A »,£12.5»2X.7HSi»^AN& - , E 1 2 . 5» 2X, 7HV w • , E 1 2 . 5# 2 X , 7HR V W AL ■♦Eia. 

l»i>/2X»7HNOITER»» 1 12, 2X» THEkKUR «, E12. & ) 

81 FORMAT (/2X>FHMISTAKE=,E12.5,37H MISTAKE-ABSOLUTE VALUE OF l.-PE/ 
IPlO/iOX, 106HY0U HAVE MADE AN ERROR IN YOUR INPUT SUCH THAT THE RAT 
210 OF STATIC TO TOTAL PRESSURE IS GREATER THAN UNITY/llX, 66HERRQR 
3C0ULD INVOLVE EITHER OR ALL OF THE F OLLOW I NG : X Ma, PT l, PE , WAV E / / 1 2 X, 
A97HTHE PROBABILITY IS VERY HIGH THAT YOUR SPECIFIED VALUE OF PE(1) 
5 IN iNAM3 IS NOT CORRECT FOR YOUR / , 1 2x , AOHS EL E TE 0 VALUES OF XMA,P 
6T1,WAVE IN $NAM2. //,12X,27HPE/P10 CANNOT EXCEED UNITY,/) 

82 F0RMaT(1H1,2X,5H?NaM1,//,2X,8HIGE0M », 112, 2X, 8HXEN0 -,E12.5,2X, 

IbHlE -,I12,2X,6HXK » , E 12 . 5 , 2 X, 8HD E T A1 El 2 , 5, 2 X, 8 HXMA 

2 £12.5, /, 2X, 8HP ri *, E 1 2 . 5 , 2 X, 6HTT1 » , E 1 2 . 5 , 2X , 8 H I GA S »,I 

312,2X,8HVIS1C1 «, £12. 5, 2X, 8HVIS1C2 « , E12 . 5, 2X, 6HV IS 2C 1 »,E12.5,/,2 
AX,8HVIS2C2 -, E12 .7,2X, BHG E 12 . 5, 2x, faHK • , El 2 . 5, 2X , 8HPR 

5 *,E12.5, //,2X,8HIctODY », 112, 2X, BHWAVE », E 12 . 5, 2X, 8H PHI I 

6-, E12 .5, 2X, 8HJ » , 1 1-: , 2 X , 8 HW », 1 12, 2X, 8H I E NT R 0 -,I12,//,2 

7X,8H5iT -, E12 . 5, 2X, 8HSMXTR -, E12 , 5, 2X, 8HK0DVIS -, II 2, 2X , 8HKTC0 
80 -,112,2X,8HTLNG1H » , E 1 2 . 5, 2 X, 8H I Y I NT - , 1 12 , / , 2X , 8HK ODAM P «,I12 

9,2X,6HXT1 », El2.5,2X,bHXT2 -, E 12. 5, 2X, 6HXT3 -, E 12 . 5, 2X, 8H 

IXTA »,E12.5,2X,8HXT5 - , E 12 . 5 , / , 2X, 8HXT6 *, E12 . 5, 2X , 8HPR T 

2 «,£12.5,2X,3HKaDPkT - , 1 1 2 , 2 X, 8HNUMtJl -, 1 12, 2X , 8HGLAR »,E12.5 
A, /, 2X, 8HPRTAR , / , ( 1 OE 12 . 5 ) ) 

83 FURMAT(/,2X,9HIENDl -, II 2, 2X, 8HPR0INC -, E12 . 5, 2X, 6HPKNT INC ■ , E12 . 5 
1,2X,8HIPR0 »,I12,2X,6HiPkNT - , 1 12 , 2x, 8HN AUX PkQ- , 1 12, / , 2X , 8HFT 

1 -,tl2,5,?X,8HK00e «, I 12,2X,8HK0UWAL -, 1 1 2 , 2 X , 8H VE L E OG -,E12,5 

1,2X,8HC0NV »,E12.5,2X,bHMTMAX », i 1 2, / , 2X, 8HK0UUN I T-, 112, 2X, 8 HIT 
2MAX *, I12,2X,9HC0NVE = , E 1 2 . 5 , / / , 2 X , 11 9HTH I S COhPLETtS THE OUTPUT 

3 OF $NAM1 wITH THE EXCEPTION uF PROVAL AND PRNTVAL. THESE VALUES A 
ARE PRINTED JUST PRIOR TO THt , / , 2 X, 22H IN IT I AL STATION PRINT.) 

6 A F0RMAT(2X,6HPRGVAL,/, (10E12.5) ) 

85 F0RMaT(2X,8HPRNTVAL , / , ( 1 OE 1 2 . 5 ) ) 

l*,/35X,lH+»10X,23HVAklBLE ENTROPY CALCULATIONS, lOX , 1H+, / , 35 X, IH* , 1 
2iX,6HlTMAX-, I 8 , 2 X , 5 HNO i T = , I 6, 8 X, 1H + , / , 3 5X, 1 H* , A 8X , 1 H*, / , 3 5X , SOH*** 

67 FQRi1aT(/,2x,63HPRINT STATIONS DEsIGNATED IN SNAMl INPUT AND GENERA 
ITEO IN SETUP, / ) 

88 FORiiAK /,2X, 125HFQr< YOUR VAKIBlE ENTROPY CASE (IENTRO-2) THE FIRST 
1 PASS OVER THE BODY IS LUOIVAlENT TO A CONVERGED CONSTANT ENTROPY! 
2IENTRD«l)/,2X,92HSGLUnON. VARIbLE ENTROPY EFFECTS ARE TREATED IN 
3SudSEuUENT PASSES; THAT IS FOR NOI T- 1 ,2, . . . ) 

69 FORMAT! /,2X, 117HYDUR VAkIBLE ENTROPY CASE HAS NOT CONVERGED FOR YO 
lUR INPUT VALUES JF CQNVE AND ITMAX. YOU MUST EITHER INCREASE ITMA 
2X, /, 2X,66H0» REDUCt YOUR CONVERGENCE LEVEL BY INCREASING THE VALUE 
3 OF CONVF . ) 

90 FORMAT! /, 2X, 77HY0UR VARIBLE CaSE ENTROPY IS CONVERGED TO YOUR SELE 
ICTED INPUT VALUE Of CQNVE.) 

91 FORMAT! 2X,A7H30UNDARY-LAYER SEPARATION INDICATED BY SOLUTION, /, 2 X, 
12HS-,E12.5,2X,5HT4UD»,E12 . 5 , 2X, 4HQSD - , E 12 . 5 ) 
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93 FG«r1Al ( 1H1» 2Xf 51HF CK J«i DiSP IS DtFINEO aS(SEE AIaA PAPER NO. 69- 
1667) FOLLOW5»//»i*X#4iH0l5P*RMi*(-l. + SQKT(i. + 2.tDISP(J*O)/RMl)))/// 
2»2X,i09HFDP YOUR CURRENT S-STaTIQN 0ISP(J«0) IS N t G AT I V E ( AL LOR E D ) 
3ANO 2.*DISP( J*0) /RhI IS less THAN -1 . (NOT ALL OrED ) ; W f 2 11 2HC ONSE 
AQUENTLY^THE VALUE PRlNrtD ABOVE IS D1SP(J»0) WHERE DISP(J*0) IS TH 
5E rWCi DIMENSIONAL uEFINlTiON. THE PR OBLEH W > 2X> 11 IHNORM ALL Y OCCURS 
6 FOR (l)EXTRFNELY FINE DS NEAR S«C. FOR BLUNT BODIES WHERE RMI IS 
7VEKY SMALL; (2) FLOW GvER VEk YW> 2 X» 1 lOHSLENDER BODY OF REVOLUTION 


8(N£6DLE). THE PRDBLEh IS EASILY SOLVED BY REFINED GRID DISTRIBUTIO 
9N. HOWEVER, IF Y OJ, /, 2X , i UOHARE HOT INTERESTED IN THE TWO ABOVE ME 
INTIJNED CASESiIGNORE THIS MESSAGE AND PROBLEM WILL CURE ITSELF,/, 2 
2X,50HAS S INCREASES h OR ALL PHYSICALLY REALISTIC FLOWS.) 

END 
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Subroutine TURBLNT 

Subroutine TURBLNT calculates either the eddy viscosity or mixing- length param 
eters required for the transitional and/or turbulent boundary- layer equations. The 
flow chart for subroutine TURBLNT is as follows: 
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The program listing for subroutine TURBLNT is as follows: 


SUBROUTINE TURBLNT ( T 3 # X L 3# F Z3» R AT03# Y> E P3# F3# EH3f VAR A, V ARB » VAR C , V 
lARDfVARE^KODAMP, lYINT, YEDGE»TAUP) 

DIMENSION VARA(JK)# VARB(JK)# VARC(JK)# VAPD(JK), VARE(JK), Y(JK)> 
1 T3(JK)/ XL3(JK), FZ3(JK), RAT03UK)# EP3(JK># F3(JK), EH3(JK), TA 
2UP( JK) 

DIMENSION PRTAR(JI), GLAR(JI), DVT(JI,1) 

EQUIVALENCE ( PRT AR# DVT < 1 # 1 ) ) 

COMMON /TRBULNT/ S» KS TR, TLNGTH# TRF AC T, D I S INC# XTl# XT2# XT6# XT3# XT A, X 
IT5#PRTW#RE#UE#XNUE# J#RMI#EPS# JP0INT#IE#WW1#WW2#WW3#WWA#WW5# NEDGE 
2#K0DVIS#A#XBE#X#PR#K0DPRT# PRT# PRTAR# GL AR#NUMB1# XK 
COMMON /UNIT/ V I SIC 1# V IS 1C 2# VI S2C1# V IS2C 2# PTl# TTl# WAVE# R# PH 10# D S# S 
1ST#RT1#P1#T1#R1#U1# AA1#TREF# VI SREF# PESTAR# TEST AR#RE STAR# UE STAR# MUE 
2STAR#YESTAR#THETA#TAUD#QSD#HD#UPLUS#DISP#PE#Z#TW#QW#RVWALD# PRO INC# 
3PRNTINC#ZS#RS 
GO TO (1#2)# KSTR+1 

1 KSTR-1 
PRTW-PRT 
STR-S 

XLAHDA-STR^(TLNGTH-1. ) /(SORTt (AL0G(50.) ) /.A12) ) 

2 CONTINUE 

CALCULATE STREAMWISE INTERMITANCY 

SN0RM-.A12*( ( ( S-STR)/XLAMDA)*^2» 

TRFACT-1. 

IF CSN0RM.LT.20. ) TRF ACT-1 ,-EXP<-SNORM) 

CALCULATE EDOY-VISCOS ITY 
2 LAYER MODEL 


IFC-0 

A2-RE*RE+UE*UE*(RMI**J )+XL3{l)*FZ3(l)/(SQPT(2,*X)*XNUE+EPS) 
A3-RE+RE*UE*UE»(RMI**J)/ { A*A*EPS*EPS*EPS*XNUE*SQRT(2.*X) > 
A4»PE+UE+XT2*DISINC/(EPS*EPS+XNUE*A) 
A5-SQPT(A2/(XL3(1)*XL3(1)*T3(1)**3)) 

NDAHP-0 
DO 8 N-2»IE 

ERA-XT3*( (YIN) / YEDGE)-XT<f) 

CALL ERF (ERA#ERB) 

YINTER-.5*(1.-ERB) 

IF (lYINT.EO.l) YINTER-1. 

IF (IFC.EQ.l) GO TO 3 

YPLUS«Y(N)^SQRT{A2/(XL3(N)*XL3{N)*T3(N)+*3)) 

IF (K0DAMP.EQ.2) YPLUS -Y ( N ) ♦ A5 
DAMP-1. 

82 — YPLUS/XT6 

IF CNDAMP.EQ.O) DAMP-1 .-EXP( B2) 

IF (DAMP. GT.0.999<J. AND. NDAMP.EO.O) NDAMP-1 

VARA(N)-A3*(RAT03{N)**J)+ABS(FZ3(N))/(XL3(N)*T3(N)*+3) 

VARB(N)-DAMP 

IF (KODVIS.EQ.l) GO TO A 
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XMIXL-XT1*A*Y(N)*DAMP»EPS 
EPl-l.*TRFACT*VARA(N) + XMIXL+Xf^IXL 

VARC (N)»EP1 
3 CONTINUE 


OUTER EDDY VISCOSITY LAW 


EP2"l.+TRFACTtAA + YINTERMXL3(N)*T3(N)*T3(N)) 
VARD(N)-EP2 
IF (IFC.EO.l) GO TO 6 

IF (EP1.LE.EP2.AND.IFC.EO.0) GO TO 5 
IFC-1 
JPOINT-N 
GO TO 6 

MIXING LENGTH MODEL 

A CONTINUE 
YRAT-1. 

IF <Y(N).LE.YEDGE) YR AT- Y ( N ) / YEDGE 
XMIXL-XT5*TANH( XT1^YRAT/XT5)*EPS+A*YEDGE*DAMP 
EP3(N)-1.+TRFACT*YINTER*VARA(N)*XMIXL*XMIXL 
GO TO 7 

5 EP3(N)-EP1 
GO TO 7 

6 EP3CN)-EP2 

7 CONTINUE 
VARE <N)-XMIXL/YEOGE 

8 CONTINUE 

CALCULATE TURBULENT PRANDTL NUMBER 

DO 12 N-2#IE 
GO TO (ll^Q^lO)# KODPRT 

9 PRT»PRTW-PRTW*.5*( (Y<N)/YEDGE)+*2) 

IF (Y(N).GE. YEDGE) PRT-PRTW+.5 
GO TO 11 

10 GL-Y<N)/YEDGE 

IF (GL.GT.l) GL-1. 

CALL lUNI ( JI,NUMB1#GLAR,1»DVT#2#GL#PRT» IPT>IERR) 
11 EH3(N)-(PRT+(EP3(N)-1. )*PR)/(PR*PRT) 

12 CONTINUE 

DO 13 N-2»IE 

IF (EP3(N).LT.l.) EP3(N)-1. 

13 CONTINUE 
VARA(1)-VARA(2 ) 

VARB(l)-0. 

VARC(1)-1. 

VAP0<1)-VARD(2) 

VARE{l)-0. 

RETURN 

END 
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Subroutine VARENT 

Subroutine VARENT reads the variable-entropy input (coordinates for shock wave; 
see fig. 2 ), computes the derivatives dRg/dz, and writes the input and derivatives. 
The flow chart for subroutine VARENT is as follows: 
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The program listing for subroutine VARENT is as follows: 


SUBROUTINE VARENT ( RR S # Z Z S# DRSDZ S, NUMBE R# KOOUNIT> A ) 

DIMENSION RRS(JL)> ZZS(JL)* DPSDZS(JL) 

NAMELIST /NAM3/ NUMBER # R P S» ZZ S» ORS DZ S 
READ (5/NAM3) 

NUMMl-NUMBER-1 
DO 1 I»2»NUMM1 

1 DRSD7S(n»(RRSn-H)-RRS( I-1))/(ZZS (I+1)-ZZS( I-in 
0RSDZSm»<RRS(2)-RRS(l) ) / (ZZS<2)-ZZS(1) ) 

OR SDZS( NUMB ER )■( RRS ( NUMB ER )-RPS( NUMB ER-1))/(ZZS( NUMBER) -ZZS( NUMBER 
1 - 1 )) 

WRITE(6#3) NUMBEP#RPS 
WRITE(6,<») ZZS 
WRITE(6#5) DRSDZS 
C5-1. 

IE (KODUNIT.EQ.I) C5-3 .28083985 
C6-C5/A 

DO 2 I«l, NUMBER 

RRS(I)-RRS(I)*C6 

ZZS(I)-ZZS(I)*C6 

2 CONTINUE 

3 F0RMAT(2X# 5H$NAM3» //> 2X, 8HNUMBER »#I12#/#2X,3HRRS#/#(10E12.5)) 

4 FORMAT(2X,3HZZS,/#aOE12.5>) 

5 F0RMAT(2X#6H0RSDZS# /# ( 10E12.5) ) 

RETURN 

END 



APPENDIX C 
Subroutine SOLVE 


ORIGINAL PAG£ m 
OF POOR QUAUTY 


Subroutine SOLVE is the main subroutine of program VGBLP wherein the nonsimilar 
laminar, transitional, and/or turbulent boundary- layer equations are solved using a 
coupled algorithm for the energy and momentum equations. The continuity equatiorv is 
using the trapezoidal rule. An iterative loop is provided to assure conver- 
gence of the system of equations to a preselected value. The flow chart for sub- 
routine SOLVE is as follows: 


Q SOLVE ^ 

i 

Set boundary conditions 



Make first guess at 


profiles 

i 

Set up matrix elements 

for momentum and energy 
equations 

y 

Calculate solution vector 


± 

Integrate 

continuity equation 

i 

Compute 

normal derivatives 


Y 

Update viscosity 
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The program listing for subroutine SOLVE is as follows: 


SUBROUTINE SOLVE ( KOD WAL # ARGM, TR, V IS2C 2» X AL, XB E# X3# DX2> TW» I GAS# IE# 
1J#NIT#RVWAID#R1#U1#WW1#WW2#WW3#WWA#WW5#XNUE# XK# EPS# M# I BODY# FT) 
COMMON /SOL VI/ FI ( JK ) # F2 < JK ) # F3 UK > # T1 ( JK ) # T2 ( JK ) , T3 ( JK ) # VI ( JK ) # V2 
1UK)#V3UK)#EP2UK)#EP3(JK)#FZ2(JK)#FZ3(JK)#TZ2(JK)#TZ3(JK)#XL2(JK 
2)#XL3( JK)»XLP2( JK)#XL P3 ( JK ) # RATOl (JK )# R AT02 ( JK )# R AT03 ( JK ) # EH2 UK ) # 
3EH3UK)#DRAT01( JK)#DRAT02UK)#DRAT03( JK)#VARA( JK)#VARBUK)#VARC ( JK 
A)#VARDUK)#VARE( JK)#YUK )#EP1( JK)#EH1( JK)#XL1( JK)#XLP1UK) 

COMMON /S0LV2/ Z 1# Z2# Z 3# Z <t# Z5# TE# F AA# F AB#F AC# F AD# BEXl# BEX2# BEX3# B E 
IX^ 

COMMON /ME SHI/ XN ( J K ) # DN ( JK ) # Y1 ( JK > , Y2 UK ) # Y3 C JK ) # YA ( JK ) # Y5 ( JK ) # Y6 
l(JK) 

DIMENSION XKl(JK)# XK2(JK)# XK3UK)# XMl(JK)# XM2UK)# XM3<JK) 

THIS SUBROUTINE SOLVES THE COMPLETE B.L. EQUATIONS 

SET UP THE BOUNDARY-CONDITIONS 
XKl(l)-XK2(l)»XK3(l)«XM2(l)-XM3(l)-0. 

IF (KODWAL .EQ.l) GO TO 1 
TZ3(1)»ARGM/XL3(1 ) 

GO TO 2 

1 CONTINUE 
XM1(1)-TW/TE 
T3(1)-XM1(1) 

2 CONTINUE 
F3(l)-0. 

IF (NIT.GT.l) GO TO 6 

MAKE FIRST GUESS AT THE CURRENT STATION PROFILES 
DO 5 N-2#IE 

T3(N)»ZAAT2(N)-Z5*T1(N) 

F3<N)-ZAAF2(N)-Z5+F1(N) 

IF <T3(N) .LT.O. » WRITE (6#18) T3(N)#N#M 
IF CT3(N).LT.O.) T3 { N ) • , 5A ( T2 ( N ) +T1 ( N ) ) 
RAT03{N)-ZA*RAT02{N)-Z5*RAT01(N) 

IF (IGAS.E0.2) GO TO 3 

XL3IN)-(<l.+TR)*SQRT(T3tN) )/(T3(N)*TR) ) 

XLP3(N)»XL3(N)A(TR-T3(N) ) /(2.*T3(N)*(T3(N)+TR) ) 

GO TO A 

3 CONTINUE 

XL3<N)-T3(N)**(VIS2C2-1. ) 

XLP3(N)-(VIS2C2-l.)A(T3(N)*^<VIS2C2-2.) ) 

A CONTINUE 

TZ3(N)-TZ2(N) 

FZ3(N)-FZ2<N) 

EP3<N)-EP2CN) 

V3<N)-2A+V2(N)-Z5AV1(N) 

EH3(N)-EH2 (N) 

5 CONTINUE 

6 CONTINUE 
IM-IE-1 
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C 

C SET UP MATRIX ELEMENTS 

C 

DO 8 N-2#IM 
C 

C MOMENTUM EQUATION 

C 

A11«X3*F3IN)/0X2*FT 
A12-V3(N)/(DN(N)4DN(N-1) ) 

CLEBP-.5*RAT03(N+1)*^(2*J)^XL3(N+1 )*EP3(N41)4.5*RAT03(N)^*(2*J) +XL 
13(N)*EP3(N) 

CLEBM».5*RAT03(N-1)**(2*J)*XL3(N-1 )*EP3(N-1)4.5*RAT03(N)*+(2*J) ♦XL 

13<N)+EP3(N) 

A13-XBE 

A14--XBE+F3(N)^+2 
Al— A12-Y3(N)^CLEBM 

Bl-All+Zl4Yl(N)+CLEBP4Y3 {N)+CL£BM42.*A13+F3(N) 

C1-A12-Y1CN)^CLEBP 

Dl-0. 

El— A13 
Hl-0. 

G1«A11+(Z2*F2(N)-Z3^F1(N) )-AlA 
C 

C ENERGY EQUATION 

C 

CLEHP»*5*RAT03(N4l)^+(2^J)^XL3(N4l)+EH3(N+l)4*5*RAT03(N)^^(2^J)*XL 

13(N)^EH3(N) 

CLEHM».5*RAT03(N-1) ♦♦(2* J )+XL3{ N-1 l+EHBIN-l )+,5+RAT03(N)+^(2^J)*XL 
13(N)+EH3<N) 

A15»XAL*XL3(N)^EP3(N)»RAT03(N)^+(2^J) 

A16»-A15^FZ3(N)%^2 
A17-FZ3(N)/(0N(N-1)4DN(N) ) 

A2-2.^A15+A17 

B2»0. 

C2 — A2 

02— A12-Y3(N)+CLEHM 
E2-A11^Z1+Y1(N)+CLEHP4Y3(N)^CLEHM 
H2-A12-Y1(N)^CLEHP 
G2-A11+(Z2+T2(N)-Z3+T1(N) )4A16 
IE (KODWAL.EQ.l ) GO TO 7 
IF (N.GT.2) GO TO 7 

DID«(C2*Dl-Cl^D2)-( (C2*Hl-CltH2)+( (<l.+XK)+^2)-l.) ) 
XM1{1)-{(C2+G1-C1^G2)4<C2*H1-C1*H2)+(XK+(1,4XK>+DN(1))+TZ3(1))/DI0 
XM2(D — (C2^B1-C1 + B2)/DI0 

XM3(1> — ( (C2*El-Cl + E2)+( (C2 + Hl-Cl + H2)^( (l.+XK ) ♦♦Z) ) ) /DID 
7 CONTINUE 

C 

C SET UP MATRIX ARRAYS 

C 

B1S-B1+A1^XK2(N-1 )+Dl+XM2(N-U 
B2S»B24A2^XK2(N-1 )4D2+XM2(N-1) 

EIS-E14A1+XK3(N-1)401^XM3CN-1» 

E2S»E2+A2^XK3(N-1)402^XM3(N-1) 

GlS«Gl-Al+XKl(N-l)-Dl+XMltN-l) 
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G2S»G2-A2*XK1{N-1)-D2*XM1(N-1) 

D-1,/(B1S*E2S-E1S + B2S ) 

XKl (N)-D* (G1S*E2S-G2S*E1S) 

XK2(N)«D>ME1S*C2-C1*E2S) 

XK3(N)-D*(E1S*H2-H1»F2S) 

XMl (N)-D*(B1S*G2S-B2S+G1S) 

XM2<N)-D+(C1+B2S-B1S+C2) 

XM3(N)-D* (H1*B2S-B1S*H2) 

8 CONTINUE 
KON-IM 

CALCULATE THE SOLUTION VECTOR 
00 9 N-2# IM 

F3(K0N)-XK1(K0N) +XK2(K0N)*F3(K0N+1 )+XK3(K0N)+T3(K0N+l) 
T3(K0N)»XM1(K0N) + XM2(K0N)*F3(K0N + 1 »-*-XM3(K0N)*T3(K0N + l) 

KON-KON-1 

9 CONTINUE 

IF (KODWAL.EQ.l) GO TO 10 

T3(1)-(XK*(1.+XK)*DN<1 ) + TZ3(l)-(l.+XK)+*2*T3(2)+T3C3) ) /n.-a. + XK) 
1**2> 

TW-T3(1)*TE 

10 CONTINUE 

INTEGRATE CONTINUITY EQUATION 

V3(1)-(RVWALD*FAA)/(R1»U1»EPS*XNUE) 

BO-0. 

A28-.5*X3/DX2#FT 
DO 11 N-2#IE 

B«A28+CZ1*F3(N)~Z2*F2(N) ♦Z3*F1(N) )*F3(N)*,5 
V3(N)«V3(N-1 )-(B+B0)*DN(N-l) 

BO-B 

11 CONTINUE 

COMPUTE NORMAL DERIVATIVES 

DO 12 N-2>IM 
RDN-1./(0N(N-1)+DN(N) ) 

TZ3(N)-(T3(N*1)-T3(N-1) )*RDN 
FZ3CN)-(F3(N*l)-F3(N-in ♦RDN 

12 CONTINUE 
FZ3CIE)-0. 

TZ3(IE)-0. 

IF (KODWAL.EQ.l) TZ3 ( 1 ) • (-WW1»T3 ( 1 ) +WW2+T3 ( 2 )-WW3^T3 ( 3 ) +WW4*T3 ( ^ ) ) 
1/(WW5*DN(1) ) 

FZ3(1)«(-WW1*F3(1)+WW2*F3(2)-WW3*F3(3)+WWA + F3(A) )/(WW5>*'DN(l ) ) 

UPDATE VISCOSITY 

KON-IE+2 
DO 17 N-1#K0N 
IF (T3(N).LT.O. ) GO TO 13 
GO TO 1<* 
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13 WRITE(6fl9) T3(N)#N#M 
STOP 500 
CONTINUE 


ORIfiif'iAL PA--^s t3> 
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IF (IGAS.E0.2) GO TO 15 
XL3(N)-( (l.+TR)*S0RT(T3(N) )/{T3(N) +TR) ) 
XLP3(N)-XL3<N)*<TR-T3<N) ) / ( 2 • *T3 (N ) ♦ ( T3 ( N ) +TP ) ) 
GO TO 16 

15 CONTINUE 

XL3(N)-T3CN)*+(VIS2C2-1. ) 
XLP3(N)-(VIS2C2-1,)*(T3(N)**(VIS2C2-2. ) ) 

16 CONTINUE 

IF (IBODY.NE.l) GO TO 17 
IF (M.NE.2) GO TO 17 
T1(N)»T3<N) 

F1(N)-F3(N) 

RAT01(N)-RAT03(N) 

DRAT0KN)-DRAT03(N) 

XL1(N)-XL3<N) 

XLP1(N)-XLP3(N) 

V1(N)-V3(N) 

17 CONTINUE 
RETURN 


18 FORMAT (2X,12HNEGATIVE T 3 F 1 5 . 7, 2 X, 27H REPLACE WITH MEAN OF T1,T2 
1,2X»7HAT N - ,I5#2X#20H AND AT STATION M • ,15) 

19 FORMAT!//, 2X,125HS0LUTI0N HAS RESULTED IN GENERATION OF NEGATIVE T 
TEMPERATURE WHICH CANNOT BE PHYSICALLY ACCEPTED. THIS GENERALLY OCC 
2URS IN THE,/,2X,125HREGI0N OF ADVERSE PRESSURE GRADIENT AND/OR HAS 
3S INJECTION AT WALL BOUNDARY AND IS AN INDICATION OF BOUNDARY LAYE 
AR SEPARATION., //, 2X, 118HIF DPEDS IS NEGATIVE AND THERE IS NO MASS 
5INJECTI0N AT THE WALL BOUNDARY, THE PROBLEM IS PROBABLY CAUSED BY 
6T00 COURSE, /,2X,32HA STEP SIZE IN THE S-C00RDINATE.,//,2X,6HT3(N>- 
7,F15.7,2HN-,I5,2X,2HM-,I5) 

END 
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_ .A '#^1 E 


Subroutine MESH generates the grid in the y-coordinate and obtains the difference 
quotient coefficients descriptive of the grid. The flow chart for subroutine MESH is 
as follows: 



The program listing for subroutine MESH is as follows: 


SUBROUTINE MESH { K> ETAMAX> IE# IGEOM, DETAl ) 

PEAL K 

COMMON /MESHl/ XN(JK)#DN(JK)#YKJK)#Y2(JK)#Y3(JK)#YA(JK)#Y5(JK)#Y6 
1( JK) 

XN(l)-0. 

JKPl-JK 

IGEOM-l SPECIFY ETAMAX#IE#K #»2 SPECIFY ETAMAX# IE# OETA ( 1 ) 

IF (K.EQ.l.) GO TO 5 
IF nGE0M,EQ.2) GO TO 1 
GO TO A 
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1 CONTINUE 
XKO-1.0001 
DIF-1. 

RATIO-ETAMAX/DETAl 

NIT-0 

2 CONTINUE 

All-1. +(XK0-1. )*RATI0 
A12-1./FL0AT(IE-1) 

XKN-A11++A12 

NIT-NIT+1 

DIF-ABS(1.-XKN/XK0) 

XKO-XKN 

IF (NIT.6T.20) OIF-0. 

IF (DIF. GT. 0.0005) GO TO 2 

NIT2-0 

DIF2-1.0 

3 CONTINUE 
ANl-FLOAT(IE-l) 

AN2- ANl-1. 

FK-(XKO)*+AN1-1.-(XKO-1. )+RATI0 

DFK-AN1+XK0++AN2-RATI0 

XKN-XKO-FK/DFK 

NIT2-NIT2+1 

DIF2-ABS(1.-XKN/XK0) 

XKO-XKN 

IF (NIT. GT. 100) DIF2-0. 

IF (DIF2. GT.O. 00000001 ) GO TO 3 
WRITE (6»9) NIT#NIT2/DIF,DIF2#XKN 
K-XKN 

A CONTINUE 

DN(1)-((1-K)/(1-K*^(IE-1 ) ) )*ETAMAX 
GO TO 6 

5 DN(1)-ETAMAX/(IE-1) 

6 XN(2)-DN(1) 

DN(2)-K*DN(1) 

DO 7 N-3»JKP1 
DN(N)-(K*+(N-1) )+DN(l) 

7 XN(N)-XN(N-1)+DN(N-1) 

DO 8 N-2,JKP1 
Dl-DN(N-l) 

D2-DN(N) 

Yl(N)-2./(D2+(Dl+D2)) 

Y2(N)-2./(Dl*D2) 

Y3(N)-2./(Dl^(01+D2)) 

YA(N)-D1/(D2*(D1+D2)) 

Y5(N)-(D1-D2)/(D1^D2) 

Y6(N)-D2/(D1"*'(D1 + D2) ) 

8 CONTINUE 
RETURN 


FORMAT (2X#2I5»3E16,8) 
END 
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Subroutine SIMILAR 

Subroutine SIMILAR generates the similar solutions for the boundary- layer equa 
tions at the initial station = 0) . The flow chart for subroutine SIMILAR is as 
follows : 
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The program listing for subroutine SIMILAR is as follows: 


ommm. PMm vs 

OF POOR QL’ALITV 


SUBROUTINE SIMILAR ( I E# I GAS # XAL » X8 E» PR# KODW AL , TW# DTDZW, XK, TR , V I S2C 
12# F3#F2#T3#T2#V3# V2#DF0Z#0TDZ#VIS# EP#ARGM#OW) 

DIMENSION F3UK)# F2<JK)# T3(JK)# T2(JK)# V3(JK)# V2UK)# OFOZ(JK) 
1 # DTDZ(JK)# VIS(JK)# EP(JK)# E(JK)# F(JK) 

COMMON /MESHl/ XN(JK)#DN(JK)#Y1(JK)#Y2(JK)#Y3(JK)#YA(JK)#Y5(JK)#Y6 
l(JK) 

C THIS SUBROUTINE GENERATES THE SELF-SIMILAR SOLUTIONS TO THE 

C BOUNDARY-LAYER EQUATIONS 

NITMAX-30 
CONV-.OOOl 

C INITIALIZE THE PROFILES 

DO 1 N-1#JK 

F3{N)-F2(N)-T3(N)-T2(N)-VIS(N)-EP(N)-1,0 
DFDZ(N)-DTDZ(N)-0.0 
V3(N)-V2IN)--XN(N) 

1 CONTINUE 
DFDZW-1. 

NIT-1 

2 CONTINUE 

IF IKODWAL .NE.l) DTDZW-OW*ARGM/VIS (1 ) 

SET UP FOR ENERGY EQUATION 

E(IE)-0. 

F(IE)-1*0 
IM-IE-1 
N-IM 

DO 3 I-2#IM 

RDN-1,/(DN(N-1)+DN(N) ) 

VISM1-.5*VIS(N-1)+.5*VIS (N) 

VISPl-.5+VIS(N+l)*.5AVIS(N) 

A--V3(N)*RDN-VISM1/PR*Y3(N) 

B-VISP1/PR*Y1 (N)^VISM1/PP+Y3(N) 

C-V3 (N)+RDN-VISP1/ PR+Yl ( N) 

D-XAL+VIS IN)*DFDZ (N)**2 
ECN)— A/(B + C*E(N + 1) ) 

F«N)-(D-CtF(N+l ) ) /(B+C*E (N+1) ) 

N-N-1 

3 CONTINUE 
T3(l)-TW 

IF (KODWAL.EO.l) GO TO A 
HEAT TRANSFER BOUNDARY CONDITION 

ANMR-XK*(1*+XK)^DN(1)+DTDZW-(1.+XK )*(1.+XK)»F(2>+E(3)*F(2)+F{3 ) 
DNMR-1.-C1.+XK)^(1.+XK)+E(2)*<1.*XK)+(1.+XK)-E(3)*E(2) 

T311 )-ANMR/0NMR 
A CONTINUE 
DO 5 N-2# IE 
T2(N)-T3(N) 

T3(N)-E(N)+T3(N-1)*F(N) 
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5 CONTINUE 
DO 7 N-1, IE 

IF (T3(N) .LT.O. ) T3(N)-1. 

IF (IGAS.EQ.2) GO TO 6 

VIS(N)-({ 1.4TR)*S0RT(T3(N))/(T3(N)+TR) ) 

GO TO 7 

6 CONTINUE 

VIS(N)-T3(N)**(VIS2C2“1. ) 

7 CONTINUE 

SET UP FOR SOLVING MOMENTUM EQUATION ORIGINAL fS 

OF POOR QUALITY 

E(IE)-0, 

FCIE)-1. 

IM-IE-1 

N-IM 

DO 8 I-2#IM 

RON-1. /(ON(N-l )+DN(N) ) 

VISMl-.5*VIS(N-l)+.5+VIS(N) 

VISP1-.5*VIS(N+1)+.5*VIS (N» 

A--V3(N)*RDN-VISM1*Y3(N) 

B-VISP1^Y1(N)+VISM1*Y3 (N)+2.*XBE^F3(N) 

C-V3(N)+R0N-VISP1*Y1(N) 

D-XBE*«F3(N)**2^T3(N) ) 

E (N)— A/( B + C*E(N + 1) ) 

F<N)-(D-C*F(N + 1) )/<B+C + E(N*D) 

N-N-1 

8 CONTINUE 
F3(l)-0, 

DO 9 N-2# IE 
F2(N)-F3(N) 

F3(N)-E(N)*F3(N-1)+F(N) 

9 CONTINUE 

IF (NIT.LT.5) GO TO 11 
DO 10 N-2#IE 
F3(N)».5^(F3(N)*F2(N) ) 

T3(N)-,5*(T3(N)tT2(N) ) 

10 CONTINUE 

11 CONTINUE 

DO 12 N-2, IM 

RON-1. /(DN(N-1H-DN(N) ) 

DTDZ(N)-<T3(N+1)“T3(N-1) )*PON 
0F0Z{N)-(F3(N+1)-F3(N-1))»RDN 

12 CONTINUE 
XK1-XK<-1. 

XK2-XK1*^2 

DTDZm-DTDZW 

IF IKODWAL.EO.I ) DTOZ ( 1 ) - ( ( 1 .-XK2) ♦TBC 1) ♦XKZ+TB < 2 )-T3 ( 3 ) ) / { XK+XKl* 
IDNID) 

DFDZ(l)-( <1.-XK2)*F3(1 )+XK2*F3(2)-F3(3) )/(XK*XKl*DN(l) ) 

0TDZnE)-0. 

DFDZ(IE)-0. 



90 


111 


o o o 


APPENDIX C 


OR!GS^JAL PAGE \S 
OF POOR QUALITY 


SOLVE THE CONTINUITY EQUATION 

V3(l)-0. 

DO 13 N-2#IE 

V3(N)«V3(N-1)-DN(N-1)%,5*(F3(N)+F3(N-1) ) 

13 CONTINUE 

0IF«ABS(1.-DFDZW/DFDZ (1) ) 

OFOZWDFDZd) 

NIT-NIT+1 

IF <NIT,6T,NITMAX) WRITE (6,15) DIF, NIT 
IF (NIT.GT.NITNAX) DIF-0. 

IF (0IF.6T.C0NV) GO TO 2 
C ASSIGN PROFILES AT THE PREVIOUS STATION 

C AND PRINT THE SIMILAR SOLUTION 

WRITE (6,18) 

WRITE (6,19) 

DO lA N-1,IE 

WRITE (6,20) XN(N),F3(N),T3(N),V3(N),DFOZ(N),DTDZ(N),VIS(N) 
F2(N)-F3(N) 

T2(N)-T3(N) 

V2(N)-V3(N) 
lA CONTINUE 

WRITE (6,17) NIT 
WRITE (6,16) 

RETURN 

C 

C 

C 

15 format (2X,27HTHE ERROR IN WALL SHEAR I S-,£l 5 .7, 5HAFTER , 15 , 1 2HI TER 
lATIONS- /) 

16 FORMAT (1H0,20X,26HINITIAL STATION PARAMETERS,/) 

17 FORMAT (9X, AOHCONVERGED SELF-SIMILAR SOLUTION R EOUIR ED, 1 5, llHI T ER A 
ITIONS.,/) 

18 FORMAT (1H1,20X,A6HSIMILAR SOLUTION REQUIRED TO INITIATE MARCHING, 
IlOH PROCEDURE, ///, 20X, lAHPROFILE VALUES,/) 

19 FORMAT (15X,3HETA,11X, AHU/UE,10X,<»HT/TE,10X,1HV,13X,2HFZ,12X,2HTZ, 
112X,2HXL,//) 

20 FORMAT (10X,7E1A.6) 

END 
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Subroutine TABLE 

Subroutine TABLE reads tabular input data for body geometry, pressure distribu 
tion, and wall-boundary conditions. These inputs are nondimensionalized and dis- 
tributed to the solution stations designated in the SS array. The flow diagram for 
subroutine TABLE is as follows: 
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The program listing for subroutine TABLE is as follows: 


OmmiM, PAGE IS 
OF POOn QUALITY 


SUBROUTINE TABLE ( I ENDl# SD#R1,U1» A,TREF,KODWAl# VISREF>KODUNlT, AWT) 
DIMENSION PE(JJ), Z(JJ), RMKJJ), TW(JJ), SUJ), RVWALD(JJ), QWIJJ 
1)# PD(JH), ZED(JH), RMIDO(JH), SSUH) 

DIMENSION DVTl(JJ,3), ANS1I3) 

EQUIVALENCE ( Z* DVTl ( 1# 1 ) ) # < RM I» DVTl ( 1# 2 ) ) # ( P E» DVTl ( 3 ) ) 

NAMELIST /NAM2/ NUMBER# L# PE# Z# R MI> TW» RVWALD# OW#S» SS 
DATA Cl/.020885^346/#C2/1.8/#C5/3,280839895/#C15/.006365880A/#C16/ 
1.0000881/# SS (2 )/0./#L/W 
DO 1 I"1#JJ 

own )-RVWALD(I )«0,0 

RMim-1. 

Twn )»AWT 

1 CONTINUE 
READ (5#NAM2) 

IF (E0F<5)) 2#3 

2 STOP 5 

3 CONTINUE 

WRITE(6#23) NUMBER#L#S 
WRITE(6#2A) Z 
WRITE(6#25) RMI 

IFIKODWAL.EQ.I ) WRITE(6#26) TW 
WRITEC6#27> QW 
WRITE(6#28) RVWALD 
WRITE(6#29) PE 
WRITE(6#30) SS 
C 
C 

IF (KODUNIT.NE.l) GO TO 7 

CONVERT iNAM2 INPUT DATA TO U.S. STANDARD UNITS 

DO A I»1#IEND1 
A SS(I)-SS(I)+C5 
DO 6 I«1#NUMBER 
S(I )«S(I)^C5 
Z(I)-Z(I)*C5 
RMK I)»RMI( I)»C5 
PEm»PE(I)*Cl 
RVWALD(I)»RVWALD( n+C15 
IF (KODWAL.NE.l ) GO TO 5 
TW(I)-TW(I)*C2 
GO TO 6 

5 QW(I)«QW(I)*C16 

6 CONTINUE 

7 DO 9 I«1#NUMBER 
PE(I)-PECI)/(R1^U1+U1) 

S(I)«S(I)/A 
RMim-RMI(I)/A 
IF (KODWAL.NE.l) GO TO 8 
TWd )-TW(I)/TREF 
GO TO 9 
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8 CONTINUE 

9 Z(I)-Z(I)/A 

DO 10 I-1»IEND1 

10 SS(I)-SSU)/A 
GO TO 1? 

ENTRY TABLEl 
READ {5,NAH2) 

IF <E0F(5)) 11»12 

11 STOP lA 

12 CONTINUE 
DS-SS(l) 

IF (SS(1).GT.DS+. 000001. OR. SS(1).LT.DS-. 000001) GO TO 13 
IF CSS(2) .GT.DS+. 000001. OR. SS(2) .LT.0S-. 000001) GO TO 13 
IF (SS(3) .GT.DS-*-. 000001. OR. SS(3).LT.DS-. 000001) GO TO 13 
GO TO lA 

13 WRITE <6,22) S S ( 1 ) , SS ( 2 ) , SS ( 3 ) » DS 
STOP 77 

lA CONTINUE 
TEMP-0. 

DO 15 I-1,IEND1 
SS(I)-TEMP+SS(I) 

15 TEMP-SS(I) 

WRITE(6,31) SS 
TEMPl-0, 

TEMP2-SS(1) 

DO 16 I-1,IEN01 

SS(I)-TEMP1 

TEMP1-TEMP2 

IF (I.LT.IENDl) TEMP2-SS ( I+l) 

IPT— 1 
S0-DS*(I-1) 

IF (SS(2) .NE.O. ) SD-SS(I) 

CALL lUNl < JJ, NUMBER, S, 3, DVTl,L,SD#ANSlf IPT, lERR ) 
ZED(I)-ANSKl) 

RMIDD(I)-ANS1(2) 

PD( I )-ANSl(3) 

16 CONTINUE 
SD-TEMPl 
IENDPl-IENDl+1 
IPT--1 

CALL I UNI ( JJ,NUMBER,S,3,DVT1,L,SD, ANS1,IPT,IERR) 
ZED(IENDP1)-ANS1(1) 

RMI00(IENDP1)-ANS1(2) 

PDCIENDPl )-ANSl (3) 

SSdENDPl )-TEMPl 

WRITE (A) TW(1),RVWALD(1),QW(1),PE(1),DS 

TW0DS-2.+DS 

DO 21 I-2,IENDP1 

IPT--1 

SD-DS*<I-1) 

IF (SS(2) .NE.O.) SD-SS(I) 

IF (KODWAL.EO.l) GO TO 17 

CALL lUNI ( JJ, NUMBER, S,1,0W,L,SD, OWD, IPT, lERP) 

GO TO 18 
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17 CONTINUE 

CALL lUNI UJ#NUMBER/S#1#TW>L#S0#TWD> IPT#IERR) 

18 CONTINUE 

CALL lUNI ( JJ,NUMBER#S#1#RVWAL0#L»S0»RVWALD0#IPT#IERR) 

IF n.EQ.IENDPl) GO TO 19 

DRDZ-(RMIDD(I*1)-RMI0D(I-1>) /(ZED( I*1)-ZE0(I-1) ) 

IF (SS(2).NE.O.) TWODS-SS(I+l)-SS( I-l) 

DPE0SD»(PD< I+U-PD( I-l) ) / TWODS 
GO TO 20 

19 IF «SS(2).NE.O.) DDS-SS( I )-SS<I-l) 

DPEDSO-(PD(I )-PD(I-l) ) /DOS 

DROZ-(RMIOD(I )-RNIOD( I-l) )/(ZED(n-ZED(I-l) ) 

20 WRITE <A) SD>PD{I)#RMIDDm»TWO»ZED(I)»DPEOSD#RVWALDD#nRDZ»OWD 

21 CONTINUE 
REWIND A 
RETURN 


22 FORMAT (1X#7HSS(1) ■» F9. A# / IX, 7HSS ( 2 ) ■, F9. A, /1X,7HSS ( 3 ) ■#F9.A,/1 
1X,7HDS -,F9.A, /IX, 53HTHESE VALUES MUST BE EQUAL FOR THE STARTIN 
2G PROCEDURE) 

23 FORMAT!/, 2X,5H$NAM2,//,2X,8HNUMBER ■, I 12, /, 2X, 8HL ■,I12,/,2X, 

11HS,/,(10E12.5)) 

2A F0RMAT(2X,1HZ,/, (10E12.5) ) 

25 FORMATC2X,3HRMI,/, (10E12.5)) 

26 FORMAT(2X,2HTW, /, (10E12. 5) ) 

27 FORHAT(2X,2HQW,/, (10E12.5) ) 

28 F0RMAT(2X,6HRVWALD, /, ( 10E12.5) ) 

29 FORMAT(2X,2HPE,/, (10F12.5) ) 

30 F0RMAT(2X,2HSS, /, (10E12.5) ) 

31 FORMAT! //,2X, 117HTHE FOLLOWING SPOINT VALUES DESIGNATE THE S-COORD 
IINATE LOCATIONS WHERE THE SOLUTIONS ARE OBTAINED DUPING THE S-MARC 
2H.,/,2X,116HY0UR PRINT STATION MUST AGREE WITH ONE OR MORE OF THE 
3SP0INT LOCATIONS) IE, YOU CAN PRINT ONLY AT SOLUTION ST ATI ONS . , /, 2 
AX,113HIF the case COMPLETES WITH A NORMAL STOP AND NO OUTPUT IS PR 
5INTED, YOUR PRINT INPUT IS IN ERROR. THE ERROR CAN BE, /, 2X, 50HN0TE 
6D BY COMPARING PROVAL AND PRNTVAL WITH SPOINT. ,///, 2 X, 6HSPDINT, /, ! 
710E12.5)) 

END 
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Subroutine INUNIT 

Subroutine INUNIT converts the International System (SI) of dimensional input 
data to the U.S. Customary System of Units for calculations in the program. The sub- 
routine then converts the output data back to the SI System of Units before output. 
The flow diagram for subroutine INUNIT is as follows: 



The program listing for subroutine INUNIT is as. follows: 


SUBROUTINE INUNIT ( PROVAL# PRNTV AL> JN ) 

COMMON /TRBULNT/ S, KSTR, TL NOTH# TRF AC T» 0 1 S I NC# X Tl, XT2# XT6» XT 3# XT A, X 
1T5,PRTW,RE#UE>XNUE, J#RMI# EPS# J POINT, IE, WW1,WW2,WW3»WWA,WW5# NED6E 
2,K00VIS,A,XBE,X,PR,K0DPRT,PRT,PRTAR,GLAR,NUMB1,XK 
COMMON /UNIT/ V IS 1C 1, V IS 1C 2, VI S2C1 , V I S2C 2, PTl, TTl, WAVE, R, PHI 0, D S, S 
1ST, RTl, PI, T1,R1,U1,AA1, TREE, VISREF,PESTAR,TE STAR,RESTAR,UE STAR, MUE 
2STAR,YESTAR,THETA,TAUD,0SD,HD,UPLUS,DISP,PE,Z,TW,QW,RVWALD,PR0INC, 
3PRNTINC,ZS,PS 

DIMENSION PROVAL(JN), PRNTVALUMJ 
PEAL MUESTAR 

DATA RTl,Pl,Tl,Rl,Ul,AAl,TREF,VISREF/8*1.0/,UPLUS/0.0/ 

DATA TAUD,S,RMI,Z,ZS, RS, YE STAR, DISP, THETA, OSD, HD/1 1*0. / 

DATA CCl / .020885 A3 A6/,CC 2 / 1.8/, CCA /5. 97995/, CC 5/ 3. 280839895 /,CC 6/ . 
10019A03196/,CC7/A7,880258/,CC8/.55555555/,CC10/.167225A78/,CCll/.3 
20A8/,CC12/515.379/,CC13/113A8.93/,CC1A/20A28.0758/ 


CONVERT $NAM1 VALUES TO U.S. STANDARD UNITS 

PR0INC»PR0INC*CC5 
PRNTINC-PRNTINC>*‘CC5 
DO 1 I-1,JN 

1 PR0VAL{I)-PR0VAL(I)ACC5 
DO 2 I-1,JM 
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2 PRNTVAL (n»PRNTVAL( I)*CC5 
A»A*CC5 
U1-U1+CC5 
AA1-AA1+CC5 
SST-SST+CC5 
PTl-PTl+CCl 
Pl-Pl^CCl 
VISREF-VISREF+CCl 
TREF-TREF^CC2 
T1-T1+CC2 
TTI-TT1*CC2 
P-R+CCA 
RT1«RT1^CC6 
R1-R1+CC6 
GO TO 5 

C 

C CONVERT SNAMl VALUES TO INTERNATIONAL STANDARD UNITS 

C 

ENTRY OUTUNIT 
PROINC-PROINC*CCll 
PRNTINC-PRNTINCACCll 
DO 3 I-1>JN 

3 PROVALm«PROVAL(I)*CCll 
DO A I-1»JM 

A PRNTVAL(I)-PRNTVAL(I)*CC11 
A-A+CCll 
SST-SSTACCll 
Ul-Ul+CCll 
AAl-AAlACCll 
PT1»PT1*CC7 
P1»P1*CC7 
VISREF-VISREF+CC7 
TREF-TREF+CC8 
T1-T1ACC8 
TT1-TT1ACC8 
R-RACCIO 
PT1»RT1*CC12 
R1-R1ACC12 
60 TO 5 
C 

C CONVERT WALL VALUES TO INTERNATIONAL STANDARD UNITS 

C 

ENTRY WALLOUT 

PESTAR»PESTAP*CC7 

MUESTAR-MUESTAR+CC7 

TAUD»TAUD*CC7 

TESTAR-TESTARACC8 

S-S+CCll 

RMI-RMI+CCll 

Z-ZACCll 

ZS-ZS+CCll 

RS-RSACCll 




97 



APPENDIX C 


ORIGirmL PAGS: ws 
OF POOR QUALITY 


UPLUS-UPLUS+CCll 
DISP-DISP^CCll 
UESTAR-UESTAR^CCll 
YESTAR»YESTAR*CC11 
THETA-THETA+CCll 
RESTAR»RESTAR^CC12 
Q$D»QSD*CC13 
HO-HO+CCIA 
5 RETURN 
END 



rule 


APPENDIX C 
Function INTEGT 


CPRIG5TO: 1 

(9F tpogft 


Function subroutine INTEGT integrates the continuity equation by the trapezoidal 
. The flow diagram for function INTEGT is as follows: 


c 


INTEGT 


) 


> 


Integrate using 
trapezoidal rule 

N 

L 


c 


Return 




The program listing for function subroutine INTEGT is as follows: 


REAL FUNCTIONINTEGT(YY,NOPTS,DX) 

DIMENSION YY(NOPTS)# DX(NOPTS) 

INTEGT-0. 

IF (NOPTS.lt. 2) GO TO 2 
00 1 N«2»N0PTS 

1 INTEGT«INTEGT + (DX(N-1) /2. )*(YY(N~1)+YY(N) ) 

2 RETURN 
END 
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Subroutine SETUP 

Subroutine SETUP determines from input where profiles and wall values are to be 
printed. The flow diagram for subroutine SETUP is as follows: 


SETUP 


s 

r 

Determine stations to be 
printed as given by 
input data 

i 

r 


RETURN ^ 


The program listing for subroutine SETUP is as follows: 


SUBROUTINE SETUP (A#B#C#J#K) 
DIMENSION B(J) 

IF (A.EQ.O) RETURN 

KPLUSe-K+2 

B(K*1)-A 

IF(KPLUS2.GT. J) RETURN 
DO 1 I-KPLUS2»J 
B(I)-B(I-1)+A 
IF (B(I).GE.C) RETURN 
1 CONTINUE 
RETURN 
END 
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RESULTS FROM TEST CASES - INPUT/OUTPUT 

The input and selected output from each of the test cases is presented in the 
following order: (1) VARDIM data; (2) input data for $NAM1, $NAM2, and $NAM3 (note 

that $NAM3 data are required only for test case number 4) ; (3) initial profile 

(5=0); (4) free-stream and reference values; (5) selected wall-print locations; 

(6) selected profile-print locations. It should be noted that (5) and (6) are 
selected locations and not all the print locations obtained for the input data. The 
selected print values allow users to verify that the software has been correctly 
implemented on their computer system. 
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